/** zx_fft.c** Implementation of Fast Fourier Transform(FFT)* and reversal Fast Fourier Transform(IFFT)** Created on: 2013-8-5* Author: monkeyzx*/#include "zx_fft.h"
#include <math.h>
#include <stdlib.h>/** Bit Reverse* === Input ===* x : complex numbers* n : nodes of FFT. @N should be power of 2, that is 2^(*)* l : count by bit of binary format, @l=CEIL{log2(n)}* === Output ===* r : results after reversed.* Note: I use a local variable @temp that result @r can be set* to @x and won't overlap.*/
static void BitReverse(complex *x, complex *r, int n, int l)
{int i = 0;int j = 0;short stk = 0;static complex *temp = 0;temp = (complex *)malloc(sizeof(complex) * n);if (!temp) {return;}for(i=0; i<n; i++) {stk = 0;j = 0;do {stk |= (i>>(j++)) & 0x01;if(j<l){stk <<= 1;}}while(j<l);if(stk < n) { /* 滿足倒位序輸出 */temp[stk] = x[i];}}/* copy @temp to @r */for (i=0; i<n; i++) {r[i] = temp[i];}free(temp);
}/** 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(complex *x, int N)
{int i,j,l,ip;static int M = 0;static int le,le2;static FFT_TYPE sR,sI,tR,tI,uR,uI;M = (int)(log(N) / log(2));/** bit reversal sorting*/BitReverse(x,x,N,M);/** 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].img * uI;tI = x[ip].real * uI + x[ip].img * uR;x[ip].real = x[i].real - tR;x[ip].img = x[i].img - tI;x[i].real += tR;x[i].img += 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(complex *x, int N)
{int k = 0;for (k=0; k<=N-1; k++) {x[k].img = -x[k].img;}fft(x, N); /* using FFT */for (k=0; k<=N-1; k++) {x[k].real = x[k].real / N;x[k].img = -x[k].img / N;}return 0;
}/** Code below is an example of using FFT and IFFT.*/
#define SAMPLE_NODES (128)
complex x[SAMPLE_NODES];
int INPUT[SAMPLE_NODES];
int OUTPUT[SAMPLE_NODES];static void MakeInput()
{int i;for ( i=0;i<SAMPLE_NODES;i++ ){x[i].real = sin(PI*2*i/SAMPLE_NODES);x[i].img = 0.0f;INPUT[i]=sin(PI*2*i/SAMPLE_NODES)*1024;}
}static void MakeOutput()
{int i;for ( i=0;i<SAMPLE_NODES;i++ ){OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;}
}void zx_fft(void)
{MakeInput();fft(x,128);MakeOutput();ifft(x,128);MakeOutput();
}
程序在TMS320C6713上實驗，主函數中調用zx_fft()函數即可。
/** 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)
#endiftypedef 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_ */ [cpp] view plaincopy print?
/** 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[i].real;x[j].imag = x[i].imag;x[i].real = tR;x[i].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[i].real - tR;x[ip].imag = x[i].imag - tI;x[i].real += tR;x[i].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的復數類型
[cpp] view plaincopy print?
typedef?struct?{??
????float?real;??
????float?imag;??
}?COMPLEX;??
typedef struct {float real;float imag;
} COMPLEX;
=====================
增加：FFT頻譜結果顯示
=====================
[plain] view plaincopy print?
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');??
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');
