傅里叶变换:
傅里叶变换是一种线性的积分变换。它的理论依据是:任何连续周期信号都可以由一组适当的正弦曲线组合而成,即使用简单的正弦、余弦函数(如sinx,Acos(ωx+θ)),可以拟合复杂函数。
使用正弦曲线的原因:在信号处理中,正弦曲线可以更简单地处理信号,且一个正弦曲线信号经过处理仍是正弦曲线,只有幅度和相位可能发生变化,但是频率和波形不变。
在信号处理中,傅里叶变换(连续)是将时域信号积分,得到频域上的信息(将一条曲线拆分成正弦曲线后,各正弦曲线的振幅,图中红色部分):
傅里叶逆变换(连续)是将频域信号积分,得到时域上的信息(将各个正弦曲线合成后的曲线,图中蓝色部分):
其中,eiwt = cos(wx)+ i * sin(wx),(欧拉公式),表示复平面上的一个点。
傅里叶变换类型:
1. 连续傅立叶变换:非周期性连续信号;
2. 傅立叶级数:周期性连续信号;
3. 离散时域傅立叶变换:非周期性离散信号;
4. 离散傅立叶变换:周期性离散信号。
离散傅立叶变换(DFT):
对于计算机来说只有离散的和有限长度的数据才能被处理,所以计算机只能处理离散傅立叶变换,而其它的变换类型只有在数学演算中才能用到。关于离散傅里叶变换的公式:
离散傅里叶变换的算法复杂度为O( N2),其中xn表示合成的离散信号,Xk表示xn在时域上的变换(这里k没有固定值,但k越大,对xn的拟合曲线就越多,拟合更精确)。以下是一个例子:
图中原始信号有4个采样点,x轴表示时间(可以看出是等距采样),y轴表示振幅。而后面两个坐标,x轴表示频率个数(时间被积分了),y轴表示振幅。
离散傅里叶变换java实现如下,输入:离散信号值(实数),输出:经过傅里叶变换得到的一组频率(实数)。(Complex是一个复数类,为了方便阅读,省略Complex的实现)
- <span style="font-size:12px;">public class DFT {
- public Complex[] dft(Complex[] x) {
- int n = x.length;
- // exp(-2i*n*PI)=cos(-2*n*PI)+i*sin(-2*n*PI)=1
- if (n == 1)
- return x;
- Complex[] result = new Complex[n];
- for (int i = 0; i < n; i++) {
- result[i] = new Complex(0, 0);
- for (int k = 0; k < n; k++) {
- //使用欧拉公式e^(-i*2pi*k/N) = cos(-2pi*k/N) + i*sin(-2pi*k/N)
- double p = -2 * k * Math.PI / n;
- Complex m = new Complex(Math.cos(p), Math.sin(p));
- result[i].plus(x[k].multiple(m));
- }
- }
- return result;
- }
- }
- </span>
计算1000个采样信号用时:
快速傅里叶变换(FFT):
由于离散傅里叶变换的算法复杂度为O( N2 ),在处理较大计算量时花费时间较长,所以JWCooley和JohnTukey在1965年的文章中提出了Cooley-TukeyFFT算法,利用XN+k = Xk这一对称性(证明略),将DFT算法的时间复杂度降低到了O(NlogN ):
可以看到,式子被分解为两个更小的DFT(每个含有N/2个元素)。在这两个更小的DFT中,可以使用XN/2+m= Xm这一性质,因为k < N, m < N/2,所以m + N/2 < N,所以XN/2+m = Xm = Xmn + e-i2π m/N * Xmm,即可以递归。又因为X1 = x0(exp(-2*n*PI)=1),所以n=1为原点。
需要注意的是,N必须是2的倍数。在递归中,如果出现N是奇数,则等式不成立,需要转到DFT进行计算。
快速傅里叶变换java实现如下,输入:离散信号值(实数),输出:经过傅里叶变换得到的一组频率(实数)。
- <span style="font-size:12px;">public class FFT {
- public Complex[] fft(Complex[] x) {
- int n = x.length;
- // 因为exp(-2i*n*PI)=1,n=1时递归原点
- if (n == 1){
- return x;
- }
- // 如果信号数为奇数,使用dft计算
- if (n % 2 != 0) {
- return dft(x);
- }
- // 提取下标为偶数的原始信号值进行递归fft计算
- Complex[] even = new Complex[n / 2];
- for (int k = 0; k < n / 2; k++) {
- even[k] = x[2 * k];
- }
- Complex[] evenValue = fft(even);
- // 提取下标为奇数的原始信号值进行fft计算
- // 节约内存
- Complex[] odd = even;
- for (int k = 0; k < n / 2; k++) {
- odd[k] = x[2 * k + 1];
- }
- Complex[] oddValue = fft(odd);
- // 偶数+奇数
- Complex[] result = new Complex[n];
- for (int k = 0; k < n / 2; k++) {
- // 使用欧拉公式e^(-i*2pi*k/N) = cos(-2pi*k/N) + i*sin(-2pi*k/N)
- double p = -2 * k * Math.PI / n;
- Complex m = new Complex(Math.cos(p), Math.sin(p));
- result[k] = evenValue[k].plus(m.multiple(oddValue[k]));
- // exp(-2*(k+n/2)*PI/n) 相当于 -exp(-2*k*PI/n),其中exp(-n*PI)=-1(欧拉公式);
- result[k + n / 2] = evenValue[k].minus(m.multiple(oddValue[k]));
- }
- return result;
- }
- public Complex[] dft(Complex[] x) {
- int n = x.length;
- // 1个信号exp(-2i*n*PI)=1
- if (n == 1)
- return x;
- Complex[] result = new Complex[n];
- for (int i = 0; i < n; i++) {
- result[i] = new Complex(0, 0);
- for (int k = 0; k < n; k++) {
- //使用欧拉公式e^(-i*2pi*k/N) = cos(-2pi*k/N) + i*sin(-2pi*k/N)
- double p = -2 * k * Math.PI / n;
- Complex m = new Complex(Math.cos(p), Math.sin(p));
- result[i].plus(x[k].multiple(m));
- }
- }
- return result;
- }
- }</span>
计算1000个采样信号用时:
与DFT算法相比,速度有明显提高。