FFT算法的完整DSP实现（2）： FFT的DSP实现

只看该作者 · 2017-11-15 10:25

FFT算法的完整DSP实现（2）： FFT的DSP实现

=================
更正更正。。。
=================
上面程序中的BitReverse函数由于使用了malloc函数，当要分配的n比较大时，在DSP上运行会出现一定的问题，因此改用伪代码中提供的倒位序方法更可靠。
修正后的完整FFT代码文件粘贴如下，在实际的DSP项目中测试通过，可直接拷贝复用。

/*
* zx_fft.h
*
* Created on: 2013-8-5
* Author: monkeyzx
*/
#ifndef _ZX_FFT_H
#define _ZX_FFT_H
#include "../Headers/Global.h"
#define TYPE_FFT_E float /* Type is the same with COMPLEX member */
#ifndef PI
#define PI (3.14159265f)
#endif
typedef COMPLEX TYPE_FFT; /* Define COMPLEX in Config.h */
extern int fft(TYPE_FFT *x, int N);
extern int ifft(TYPE_FFT *x, int N);
#endif /* ZX_FFT_H_ */

复制代码

/*
* zx_fft.c
*
* Implementation of Fast Fourier Transform(FFT)
* and reversal Fast Fourier Transform(IFFT)
*
* Created on: 2013-8-5
* Author: monkeyzx
*
* TEST OK 2014.01.14
* == 2014.01.14
* Replace @BitReverse(x,x,N,M) by refrence to
* <The Scientist and Engineer's Guide to Digital Signal Processing>
*/
#include "zx_fft.h"
/*
* FFT Algorithm
* === Inputs ===
* x : complex numbers
* N : nodes of FFT. @N should be power of 2, that is 2^(*)
* === Output ===
* the @x contains the result of FFT algorithm, so the original data
* in @x is destroyed, please store them before using FFT.
*/
int fft(TYPE_FFT *x, int N)
{
int i,j,l,k,ip;
static int M = 0;
static int le,le2;
static TYPE_FFT_E sR,sI,tR,tI,uR,uI;
M = (int)(log(N) / log(2));
/*
* bit reversal sorting
*/
l = N / 2;
j = l;
//BitReverse(x,x,N,M);
for (i=1; i<=N-2; i++) {
if (i < j) {
tR = x[j].real;
tI = x[j].imag;
x[j].real = x.real;
x[j].imag = x.imag;
x.real = tR;
x.imag = tI;
}
k = l;
while (k <= j) {
j = j - k;
k = k / 2;
}
j = j + k;
}
/*
* For Loops
*/
for (l=1; l<=M; l++) { /* loop for ceil{log2(N)} */
le = (int)pow(2,l);
le2 = (int)(le / 2);
uR = 1;
uI = 0;
sR = cos(PI / le2);
sI = -sin(PI / le2);
for (j=1; j<=le2; j++) { /* loop for each sub DFT */
//jm1 = j - 1;
for (i=j-1; i<=N-1; i+=le) { /* loop for each butterfly */
ip = i + le2;
tR = x[ip].real * uR - x[ip].imag * uI;
tI = x[ip].real * uI + x[ip].imag * uR;
x[ip].real = x.real - tR;
x[ip].imag = x.imag - tI;
x.real += tR;
x.imag += tI;
} /* Next i */
tR = uR;
uR = tR * sR - uI * sI;
uI = tR * sI + uI *sR;
} /* Next j */
} /* Next l */
return 0;
}
/*
* Inverse FFT Algorithm
* === Inputs ===
* x : complex numbers
* N : nodes of FFT. @N should be power of 2, that is 2^(*)
* === Output ===
* the @x contains the result of FFT algorithm, so the original data
* in @x is destroyed, please store them before using FFT.
*/
int ifft(TYPE_FFT *x, int N)
{
int k = 0;
for (k=0; k<=N-1; k++) {
x[k].imag = -x[k].imag;
}
fft(x, N); /* using FFT */
for (k=0; k<=N-1; k++) {
x[k].real = x[k].real / N;
x[k].imag = -x[k].imag / N;
}
return 0;
}

复制代码

   另外，可能还需要您在其它头文件中定义COMPLEX的复数类型

typedef struct {

float real;

      float imag;

} COMPLEX;

复制代码

=====================
增加：FFT频谱结果显示
=====================

clc;

clear;

% Read a real audio signal

[fname,pname]=uigetfile(...

{'*.wav';'*.*'},...

'Input wav File');

[x fs nbits] = wavread(fullfile(pname,fname));

% Window

% N = 1024;

N = size(x,1);

x = x(1:N, 1);

% FFT

y = fft(x);

% 频率对称，取N/2

y = y(1:N/2);

% FFT频率与物理频率转换

x1 = 1:N/2;

x2 = (1:N/2)*fs/(N/2);  % /1000表示KHz

log_x2 = log2(x2);

% plot

figure,

subplot(2,2,1);plot(x);

xlabel('Time/N');ylabel('Mag');grid on

title('原始信号');

subplot(2,2,2);plot(x1, abs(y));

xlabel('Freq/N');ylabel('Mag');grid on

title('FFT信号/横坐标为点数');

subplot(2,2,3);plot(x2,abs(y));

xlabel('Freq/Hz');ylabel('Mag');grid on

title('FFT信号/横坐标为物理频率');

subplot(2,2,4);plot(log_x2,abs(y));

xlabel('Freq/log2(Hz)');ylabel('Mag');grid on

title('FFT信号/横坐标为物理频率取log');

% 更常见是将幅值取log

y = log2(abs(y));

figure,

plot(x2,y);

xlabel('Freq/Hz');ylabel('Mag/log2');grid on

title('幅值取log2');

复制代码