?

?

?

?

?

?

?

/****************************************************/

/****************************************************/

/****************************************************/

下面是具體說明

1、FFT：頻譜關(guān)于中間位置對(duì)稱，只需要觀察 ?0:1：N/2（這N/2+1個(gè)點(diǎn)）（時(shí)域采集N個(gè)點(diǎn)，頻域只需要觀察N/2+1個(gè)點(diǎn)）

2、MATLAB中FFT的頻譜，應(yīng)該看幅值

3、X軸頻率點(diǎn)的設(shè)置：采樣頻率為Fs,頻譜圖顯示的最高頻率為Fs/2（采樣定理）

? ：X軸頻率點(diǎn)：（0:1：N/2）*Fs/N

4、復(fù)數(shù)幅值修正

?5、

?

/****************************************************/

/****************************************************/

/****************************************************/

栗子及實(shí)踐部分

?

一、信號(hào)

%% FFT clear;clc;close all Fs=1000; % 采集頻率 T=1/Fs; % 采集時(shí)間間隔 N=2000; % 采集信號(hào)的長(zhǎng)度--采樣點(diǎn)數(shù) f1=33; % 第一個(gè)余弦信號(hào)的頻率 f2=200; % 第二個(gè)余弦信號(hào)的頻率t=(0:1:N-1)*T; % 定義整個(gè)采集時(shí)間點(diǎn) t=t'; % 轉(zhuǎn)置成列向量y=1.2+2.7*cos(2*pi*f1*t+pi/4)+5*cos(2*pi*f2*t+pi/6); % 時(shí)域信號(hào)

?

二、繪制時(shí)域信號(hào)

%% 繪制時(shí)域信號(hào) figure plot(t,y) xlabel('時(shí)間') ylabel('信號(hào)值') title('時(shí)域信號(hào)')

　　

　　

?

三、FFT變換、并繪制-幅值、實(shí)部、虛部

%% fft變換 Y=fft(y); % Y為fft變換的結(jié)果，為復(fù)數(shù)向量 A=abs(Y); % 復(fù)數(shù)的幅值（模） RE=real(Y); % 復(fù)數(shù)的實(shí)部 IM=imag(Y); % 復(fù)數(shù)的虛部%% 繪制fft變換結(jié)果（幅值，實(shí)部，虛部） figure subplot(3,1,1) plot(0:1:N-1,A) xlabel('序號(hào) 0 ~ N-1') ylabel('幅值') grid on%% 頻域只讀取0:1:N/2 subplot(3,1,2) plot(0:1:N-1,RE) xlabel('序號(hào) 0 ~ N-1') ylabel('實(shí)部') grid onsubplot(3,1,3) plot(0:1:N-1,IM) xlabel('序號(hào) 0 ~ N-1') ylabel('虛部') grid on

　　

　　

可以看出頻域中的點(diǎn)關(guān)于（N/2）對(duì)稱，所以只需要觀察（0:1：N/2）

?

四、更改相位

?

?

幅值不受影響，但實(shí)部或虛部的值，會(huì)出現(xiàn)0的情況==>看MATLAB中FFT的頻譜，應(yīng)該看幅值

?繪制半譜圖（幅值的）后--我們發(fā)現(xiàn)-幅值-相位-頻率---均和時(shí)域?qū)?yīng)不上。

==>進(jìn)行幅值-修正

五、進(jìn)行幅值-修正--并繪制圖形

%% fft變換 Y=fft(y); % Y為fft變換結(jié)果，復(fù)數(shù)向量 Y=Y(1:N/2+1); % 只看變換結(jié)果的一半即可 A=abs(Y); % 復(fù)數(shù)的幅值（模） f=(0:1:N/2)*Fs/N; % 生成頻率范圍 f=f'; % 轉(zhuǎn)置成列向量%% 幅值修正 A_adj=zeros(N/2+1,1); A_adj(1)=A(1)/N; % 頻率為0的位置 A_adj(end)=A(end)/N; % 頻率為Fs/2的位置 A_adj(2:end-1)=2*A(2:end-1)/N;%% 繪制頻率幅值圖 figure subplot(2,1,1) plot(f,A_adj) xlabel('頻率 (Hz)') ylabel('幅值 (修正后)') title('FFT變換幅值圖') grid on%% 繪制頻譜相位圖 subplot(2,1,2) phase_angle=angle(Y); % angle函數(shù)的返回結(jié)果為弧度 phase_angle=rad2deg(phase_angle); plot(f,phase_angle) xlabel('頻率 (Hz)') ylabel('相位角 (degree)') title('FFT變換相位圖') grid on

　　

放大后可以看到，此時(shí)，幅值-頻率都和時(shí)域一致

此時(shí)FFT的相位圖是雜亂無章的--不用擔(dān)心，沒有頻率處的相位是無意義的--我們只需要放大看各個(gè)（實(shí)際存在的）頻率點(diǎn)的相位即可

可以看到--f1=33Hz處為45度，即pi/4--是正確的

?

?

六、實(shí)際操作：請(qǐng)分析一個(gè)未知的采集信號(hào) (example.mat)，并確定該采集信號(hào)的頻率成分。其中， 信號(hào)的采集頻率 Fs = 2500 Hz

代碼

clear;clc;close all load('example') Fs=2500; % 采集頻率 T=1/Fs; % 采集時(shí)間間隔 N=length(y); % 采集信號(hào)的長(zhǎng)度t=(0:1:N-1)*T; % 定義整個(gè)采集時(shí)間點(diǎn) t=t'; % 轉(zhuǎn)置成列向量% 繪制時(shí)域信號(hào) figure plot(t,y) xlabel('時(shí)間') ylabel('信號(hào)值') title('時(shí)域信號(hào)')% fft變換 Y=fft(y); % Y為fft變換結(jié)果，復(fù)數(shù)向量 Y=Y(1:N/2+1); % 只看變換結(jié)果的一半即可 A=abs(Y); % 復(fù)數(shù)的幅值（模） f=(0:1:N/2)*Fs/N; % 生成頻率范圍 f=f'; % 轉(zhuǎn)置成列向量% 幅值修正 A_adj=zeros(N/2+1,1); A_adj(1)=A(1)/N; % 頻率為0的位置 A_adj(end)=A(end)/N; % 頻率為Fs/2的位置 A_adj(2:end-1)=2*A(2:end-1)/N;% 繪制頻率幅值圖 figure subplot(2,1,1) plot(f,A_adj) xlabel('頻率 (Hz)') ylabel('幅值 (修正后)') title('FFT變換幅值圖') grid on% 繪制頻譜相位圖 subplot(2,1,2) phase_angle=angle(Y); % angle函數(shù)的返回結(jié)果為弧度 phase_angle=rad2deg(phase_angle); plot(f,phase_angle) xlabel('頻率 (Hz)') ylabel('相位角 (degree)') title('FFT變換相位圖') grid on

　　

?

?

?

?

全部代碼下載地址

?

?

?

轉(zhuǎn)載于:https://www.cnblogs.com/WHaoL/p/6595132.html

總結(jié)

以上是生活随笔為你收集整理的FFT-Matlab初步实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。