c++ 傅里叶变换 分析

[shadowyl]

设原始信号f[n]由一个模拟频率为50Hz的正弦信号，以12800Hz的频率采样获得 f[i] = 311 * sin(2 * pi * 50 * i / 12800)的一个周期（N=256）的数据，取整数值（产生谐波），对信号进行谐波分析。
傅里叶变换程序如下：

// Pinyuchouqu.h
#include <complex>
#include <math.h>

#define pi 3.14159265
#define N 256
typedef std::complex<double> complex;

/*----频域抽取的FFT算法----*/
void fft(complex f[]);

/*----频域抽取的IFFT算法----*/
void ifft(complex f[]);

/*----显示信号数据----*/
void display(complex f[]);

#include "stdafx.h"
#include "Pinyuchouqu.h"

/*----频域抽取的FFT算法----*/
void fft(complex f[])
{
	int i, j, k, m, n, l, r, M;
	int la, lb, lc;
	complex t;
	/*----计算分解的级数M=log2(N)----*/
	for (i = N, M = 1; (i = i / 2) != 1; M++); 

	/*----FFT算法----*/
	for (m = 1; m <= M; m++)
	{
		la = pow(2.0, M + 1 - m); //la=2^m代表第m级每个分组所含节点数	
		lb = la / 2;              //lb代表第m级每个分组所含碟形单元数
		//同时它也表示每个碟形单元上下节点之间的距离
		/*----碟形运算----*/
		for (l = 1; l <= lb; l++)
		{
			r = (l - 1) * pow(2.0, m - 1);
			for (n = l - 1; n < N - 1; n = n + la) //遍历每个分组，分组总数为N/la
			{
				lc = n + lb;       //n,lc分别代表一个碟形单元的上、下节点编号     
				t = f[n] + f[lc];
				f[lc] = (f[n] - f[lc]) * complex(cos(2 * pi * r / N), -sin(2 * pi * r / N));
				f[n] = t;
			}
		}
	}

	/*----按照倒位序重新排列变换后信号----*/
	for (i = 1, j = N / 2; i <= N - 2; i++)
	{
		if (i < j)
		{
			t = f[j];
			f[j] = f[i];
			f[i] = t;
		}
		k = N / 2;
		while (k <= j)
		{
			j = j - k;
			k = k / 2;
		}
		j = j + k;
	}
}

/*----频域抽取的IFFT算法----*/
void ifft(complex f[])
{
	int i, j, k, m, n, l, r, M;
	int la, lb, lc;
	complex t;
	/*----计算分解的级数M=log2(N)----*/
	for (i = N, M = 1; (i = i / 2) != 1; M++); 

	/*----按照倒位序重新排列原信号----*/
	for (i = 1, j = N / 2; i <= N - 2; i++)
	{
		if (i < j)
		{
			t = f[j];
			f[j] = f[i];
			f[i] = t;
		}
		k = N / 2;
		while (k <= j)
		{
			j = j - k;
			k = k / 2;
		}
		j = j + k;
	}

	/*----将信号乘以1/N----*/
	for (i = 0; i < N; i++)
	{
		f[i] = f[i] / complex(N, 0.0);
	}

	/*----IFFT算法----*/
	for (m = 1; m <= M; m++)
	{
		la = pow(2.0, m); //la=2^m代表第m级每个分组所含节点数		
		lb = la / 2;    //lb代表第m级每个分组所含碟形单元数
		//同时它也表示每个碟形单元上下节点之间的距离
		/*----碟形运算----*/
		for (l = 1; l <= lb; l++)
		{
			r = (l - 1) * pow(2.0, M - m);
			for (n = l - 1; n < N - 1; n = n + la)   //遍历每个分组，分组总数为N/la
			{
				lc = n + lb;                  //n,lc分别代表一个碟形单元的上、下节点编号     
				t = f[lc] * complex(cos(2 * pi * r / N), sin(2 * pi * r / N));
				f[lc] = f[n] - t;
				f[n] = f[n] + t;
			}
		}
	}
}

/*----显示信号数据----*/
void display(complex f[])
{
	int i;
	for (i = 0; i < N; i++)
	{
		std::cout.width(9);
		std::cout.setf(std::ios::fixed);
		std::cout.precision(4);
		std::cout << f[i].real();
		std::cout.flags(std::ios::showpos);
		std::cout.width(9);
		std::cout.setf(std::ios::fixed);
		std::cout.precision(4);
		std::cout << f[i].imag() << 'i' << '\t';
		std::cout.unsetf(std::ios::showpos);
		if ((i + 1) % 3 == 0)
		{
			std::cout << std::endl;
		}
	}
}

// main.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"

/*----频域抽取主函数----*/
#include "Pinyuchouqu.h"
int main()
{
	int i;
	complex f[N];
	for (i = 0; i < N; i++)
	{
		// 设原始信号f[n]由一个模拟频率为50Hz的正弦信号，以128000Hz的频率采样获得
		f[i] = complex(static_cast<int>(311 * sin(2 * pi * 50 * i / 12800)), 0.0);
	}
	std::cout << std::endl << "原信号：" << std::endl;
	display(f);
	fft(f);
	std::cout << std::endl << std::endl << "FFT变换后的信号：" << std::endl;
	display(f);
	ifft(f);
	std::cout << std::endl << std::endl << "IFFT变换后的信号：" << std::endl;
	display(f);
}

运行结果如下：

原信号：
0.0000 +0.0000i 7.0000 +0.0000i 15.0000 +0.0000i
22.0000 +0.0000i 30.0000 +0.0000i 38.0000 +0.0000i
45.0000 +0.0000i 53.0000 +0.0000i 60.0000 +0.0000i
68.0000 +0.0000i 75.0000 +0.0000i 82.0000 +0.0000i
90.0000 +0.0000i 97.0000 +0.0000i 104.0000 +0.0000i
111.0000 +0.0000i 119.0000 +0.0000i 126.0000 +0.0000i
132.0000 +0.0000i 139.0000 +0.0000i 146.0000 +0.0000i
......

FFT变换后的信号：
0.0000 +0.0000i 0.0001-39728.5682i 0.0000 +0.0000i
0.0000 +27.3071i 0.0000 +0.0000i -0.0000 +18.3607i
0.0000 +0.0000i 0.0000 +12.7869i 0.0000 +0.0000i
-0.0000 +8.6700i 0.0000 +0.0000i 0.0000 -0.1501i
0.0000 +0.0000i -0.0000 +2.1999i 0.0000 +0.0000i
0.0000 +7.1615i 0.0000 +0.0000i -0.0000 +4.5055i
0.0000 +0.0000i 0.0000 +4.7312i 0.0000 +0.0000i
-0.0000 +4.7006i 0.0000 +0.0000i 0.0000 -6.5455i
0.0000 +0.0000i -0.0000 +12.5515i 0.0000 +0.0000i
0.0000 +15.0080i 0.0000 +0.0000i -0.0000 +1.3077i
0.0000 +0.0000i 0.0000 +4.6457i 0.0000 +0.0000i
0.0000 -0.7869i 0.0000 +0.0000i 0.0000 +4.3792i
0.0000 +0.0000i -0.0000 +1.4976i 0.0000 +0.0000i
0.0000 -1.0110i 0.0000 +0.0000i -0.0000 +4.3645i
0.0000 +0.0000i 0.0000 +3.5049i 0.0000 +0.0000i
-0.0000 +9.8708i 0.0000 +0.0000i 0.0000 -3.8911i
0.0000 +0.0000i -0.0000 +0.4026i 0.0000 +0.0000i
0.0000 -6.2052i 0.0000 +0.0000i -0.0000 +0.6435i
0.0000 +0.0000i 0.0000 -1.1827i 0.0000 +0.0000i
-0.0000 +7.5304i 0.0000 +0.0000i 0.0000 +11.6675i
0.0000 +0.0000i -0.0000 -3.3411i 0.0000 +0.0000i
0.0000 -0.9460i 0.0000 +0.0000i 0.0000 -18.9753i
0.0000 +0.0000i 0.0000 -5.2774i 0.0000 +0.0000i
-0.0000 -1.4888i 0.0000 +0.0000i 0.0000 +3.1466i
0.0000 +0.0000i -0.0000 +1.1237i 0.0000 +0.0000i
0.0000 +11.4144i 0.0000 +0.0000i -0.0000 +4.0522i
0.0000 +0.0000i 0.0000 -0.0139i 0.0000 +0.0000i
0.0000 -2.8643i 0.0000 +0.0000i 0.0000 +3.2044i
0.0000 +0.0000i 0.0000 -7.5589i 0.0000 +0.0000i
0.0000 +7.7567i 0.0000 +0.0000i -0.0000 +3.9821i
0.0000 +0.0000i 0.0000 +1.9450i 0.0000 +0.0000i
0.0000 -3.7617i 0.0000 +0.0000i 0.0000 +5.1762i
0.0000 +0.0000i -0.0000 +7.1931i 0.0000 +0.0000i
0.0000 +13.1961i 0.0000 +0.0000i 0.0000 -5.4580i
0.0000 +0.0000i 0.0000 +14.1624i 0.0000 +0.0000i
-0.0000 +7.4716i 0.0000 +0.0000i 0.0000 +3.4785i
0.0000 +0.0000i -0.0000 +8.0713i 0.0000 +0.0000i
0.0000 -4.7995i 0.0000 +0.0000i 0.0000 -7.6498i
0.0000 +0.0000i 0.0000 -0.0500i 0.0000 +0.0000i
0.0000 -8.5007i 0.0000 +0.0000i 0.0000 +0.2426i
0.0000 +0.0000i -0.0000 +5.1901i 0.0000 +0.0000i
0.0000 +3.6563i 0.0000 +0.0000i 0.0000 -11.1675i
0.0000 +0.0000i 0.0001 -6.9313i 0.0000 +0.0000i
-0.0000 +6.9313i 0.0000 +0.0000i 0.0000 +11.1675i
0.0000 +0.0000i 0.0000 -3.6563i 0.0000 +0.0000i
0.0000 -5.1901i 0.0000 +0.0000i 0.0000 -0.2427i
0.0000 +0.0000i 0.0000 +8.5007i 0.0000 +0.0000i
-0.0000 +0.0500i 0.0000 +0.0000i 0.0000 +7.6498i
0.0000 +0.0000i -0.0000 +4.7995i 0.0000 +0.0000i
0.0000 -8.0713i 0.0000 +0.0000i 0.0000 -3.4785i
0.0000 +0.0000i 0.0000 -7.4716i 0.0000 +0.0000i
0.0000 -14.1624i 0.0000 +0.0000i 0.0000 +5.4580i
0.0000 +0.0000i 0.0000 -13.1961i 0.0000 +0.0000i
0.0000 -7.1931i 0.0000 +0.0000i 0.0000 -5.1762i
0.0000 +0.0000i 0.0000 +3.7617i 0.0000 +0.0000i
-0.0000 -1.9450i 0.0000 +0.0000i 0.0000 -3.9821i
0.0000 +0.0000i 0.0000 -7.7567i 0.0000 +0.0000i
0.0000 +7.5589i 0.0000 +0.0000i 0.0000 -3.2044i
0.0000 +0.0000i 0.0000 +2.8643i 0.0000 +0.0000i
-0.0000 +0.0139i 0.0000 +0.0000i 0.0000 -4.0522i
0.0000 +0.0000i 0.0000 -11.4144i 0.0000 +0.0000i
0.0000 -1.1237i 0.0000 +0.0000i 0.0000 -3.1466i
0.0000 +0.0000i 0.0000 +1.4888i 0.0000 +0.0000i
-0.0000 +5.2774i 0.0000 +0.0000i 0.0000 +18.9753i
0.0000 +0.0000i -0.0000 +0.9460i 0.0000 +0.0000i
0.0000 +3.3411i 0.0000 +0.0000i 0.0000 -11.6675i
0.0000 +0.0000i 0.0000 -7.5304i 0.0000 +0.0000i
-0.0000 +1.1827i 0.0000 +0.0000i 0.0000 -0.6435i
0.0000 +0.0000i -0.0000 +6.2052i 0.0000 +0.0000i
0.0000 -0.4026i 0.0000 +0.0000i -0.0000 +3.8910i
0.0000 +0.0000i 0.0000 -9.8708i 0.0000 +0.0000i
0.0000 -3.5049i 0.0000 +0.0000i 0.0000 -4.3645i
0.0000 +0.0000i -0.0000 +1.0110i 0.0000 +0.0000i
0.0000 -1.4976i 0.0000 +0.0000i 0.0000 -4.3792i
0.0000 +0.0000i 0.0000 +0.7870i 0.0000 +0.0000i
0.0000 -4.6457i 0.0000 +0.0000i 0.0000 -1.3077i
0.0000 +0.0000i 0.0000 -15.0080i 0.0000 +0.0000i
0.0000 -12.5515i 0.0000 +0.0000i -0.0000 +6.5454i
0.0000 +0.0000i 0.0000 -4.7006i 0.0000 +0.0000i
0.0000 -4.7312i 0.0000 +0.0000i 0.0000 -4.5055i
0.0000 +0.0000i 0.0000 -7.1615i 0.0000 +0.0000i
0.0000 -2.1999i 0.0000 +0.0000i -0.0000 +0.1501i
0.0000 +0.0000i 0.0000 -8.6700i 0.0000 +0.0000i
0.0000 -12.7870i 0.0000 +0.0000i 0.0000 -18.3606i
0.0000 +0.0000i 0.0000 -27.3071i 0.0000 +0.0000i
-0.0004+39728.5682i
请教大家：变换后的实部虚部各代表什么，各次谐波如何分析（幅值和相角怎么算）？

[瓶盒]

引用 6 楼 worldy 的回复:

[quote=引用 4 楼 bottlebox 的回复:] 实部和虚部的平方和对应该点的振幅，具体如何对应参看：http://bbs.csdn.net/topics/360187658

这个说法有误，平方和再开方[/quote] 对应就可，开方增加运算量，意义不大

[worldy]

每个点代表的归一化的频率点，每个输出值就是对应频率的频域向量，模代表幅度，相角代表时延 向量的的相角应该会算吧？

[worldy]

FFT变换后的信号： 0.0000 +0.0000i 0.0001-39728.5682i 0.0000 +0.0000i 看起来红色的数据应该不对！

[Adol1111]

引用 2 楼 worldy 的回复:

FFT变换后的信号： 0.0000 +0.0000i 0.0001-39728.5682i 0.0000 +0.0000i 看起来红色的数据应该不对！

这么长一堆数据，你竟然还能挑错....

[瓶盒]

实部和虚部的平方和对应该点的振幅，具体如何对应参看：http://bbs.csdn.net/topics/360187658

[worldy]

引用 3 楼 Adol1111 的回复:

[quote=引用 2 楼 worldy 的回复:] FFT变换后的信号： 0.0000 +0.0000i 0.0001-39728.5682i 0.0000 +0.0000i 看起来红色的数据应该不对！

这么长一堆数据，你竟然还能挑错....[/quote] 因为这个数据不合理，看4L给出的理由

[worldy]

引用 4 楼 bottlebox 的回复:

实部和虚部的平方和对应该点的振幅，具体如何对应参看：http://bbs.csdn.net/topics/360187658

这个说法有误，平方和再开方