快速傅里葉變換C語言實作 模擬采樣進行頻譜分析
FFT是DFT的快速演算法用于分析確定信號(時間連續可積信號、不一定是周期信號)的頻率(或相位、此處不研究相位)成分，且傅里葉變換對應的
ω
\omega
ω是連續的，從而達到分析信號成分的目的，
具體理論參考FFT百度百科原理，
下面給出代碼分析以及模擬采樣頻譜分析的結果，
對于復信號FFT：
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
# define PI 3.14159265358979323846264338327950288 //根據運算精度或要求自行截取位數
#define q 8
#define N (1<<q) /*2的q次方點數 */
#define Fs 2560 //采樣率 Fs > 2fh
typedef float real;
typedef struct{
real Re;
real Im;
}complex;
static void Asm_Mag(complex *x,int n)
{
int i;
real Mag;
for(i=0; i<n/2; i++)
{
Mag = sqrt(x[i].Re*x[i].Re + x[i].Im*x[i].Im)/N;
printf("第%d個點：%.3f hz 幅度：%.3f \n",i,(float)i*Fs/N,Mag);
}
}
void fft( complex *v, int n, complex *tmp )
{
if(n>1) {
int k,m; complex z, w, *vo, *ve;
ve = tmp; vo = tmp+n/2;
for(k=0; k<n/2; k++) {
ve[k] = v[2*k];
vo[k] = v[2*k+1];
}
fft( ve, n/2, v ); /* FFT 偶序列 */
fft( vo, n/2, v ); /* FFT 奇序列*/
for(m=0; m<n/2; m++) {
w.Re = cos(2*PI*m/(double)n);
w.Im = -sin(2*PI*m/(double)n);
z.Re = w.Re*vo[m].Re - w.Im*vo[m].Im;
z.Im = w.Re*vo[m].Im + w.Im*vo[m].Re;
v[ m ].Re = ve[m].Re + z.Re;
v[ m ].Im = ve[m].Im + z.Im;
v[m+n/2].Re = ve[m].Re - z.Re;
v[m+n/2].Im = ve[m].Im - z.Im;
}
}
return;
}
int main(void)
{
complex v[N],scratch[N];
int k;
/* 模擬數字采樣 */
for(k=0; k<N; k++) { //F0= Fs/N、f = w/2*pi、
v[k].Re = 9999 + 570*cos(690*2*PI*k/(double)Fs);
v[k].Im = 570*sin(690*2*PI*k/(double)Fs); /* 復數信號分為實部虛部 但一般采樣為純實 虛部為0*/
}
/* FFT: */
fft( v, N, scratch );
Asm_Mag(v,N); //幅度測量
}
Fs:
其中Fs是采樣率，根據奈奎斯特采樣定律，當采樣頻率Fs大于信號中最高頻率fs的2倍時(Fs>2fh)，采樣之后的數字信號完整地保留了原始信號中的資訊，此時進行快速傅里葉變換即可得到信號的全部頻率資訊，
N:
做快速傅里葉變換FFT的點數必須是2的冪次方，
Fh:
對連續時間信號進行采樣程序中，cos/sin( ω \omega ωt)其中t被替換為kTs,即t=kTs.Ts = 1 /Fs;所以代碼中的被測信號頻率再由 ω \omega ω=2*Pi * f可得 f = ω \omega ω / 2 *pi; 此時代碼中的被測信號頻率為690hz，且幅度為570，而且包含著 9999的直流分量(0hz)，
Fo:
Fo為頻譜解析度，Fo = Fs / N ;此代碼例子中為 2560 /256 = 10 hz;這個指標需根據現實題目要求進行調整，且涉及了Fs和N 這三個量可以靈活計算組合，
complex：
FFT目標為復信號，虛部在實際采樣或者模擬采樣時虛部可以為0，且一般為0，
此FFT是時域抽選方法：DIT2 F先在時域先進行奇歐倒序，頻域輸出為正序，
結果分析：
結果明顯在0hz直流分量有幅度為9999
在690hz有幅度570的分量:
在其余頻率分量處并無分幅度，得到了被測信號的全部頻率幅度資訊，
由于FFT的頻譜結果是關于奈奎斯特頻率對稱的，所以只計算一半的點即可，
其他被測信號可自行模擬測驗驗證，
對于純實信號FFT：
此時真實幅度、代碼計算應該修改為：
sqrt(x[i].Re*x[i].Re + x[i].Im*x[i].Im)*2/N;(改為純實數時,需改為2倍)
if(i==0)
Mag = Mag /2; //直流分量不需要乘以2，所以按照乘以2計算后，要除以2
被測信號：
v[k].Re = 19999 + 570*cos(690*2*PI*k/(double)Fs) + 670*cos(590*2*PI*k/(double)Fs);.//直流分量19999 690hz分量570 590hz分量670
v[k].Im = 0; //虛部為0
為什么用cos可以思考隨機信號分析內容、sin 、cos進行FFT結果一樣
純實信號FFT結果：
可見直流分量0hz幅度為19999
在590hz和690hz均有幅度670 和570
其余頻率分量并無幅度，但在真實采樣中全頻率具有高斯白噪聲，
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/271653.html
標籤：其他