/****************************************************/
/****************************************************/
/****************************************************/
下面是具體說明
1、FFT:頻譜關於中間位置對稱,只需要觀察 0:1:N/2(這N/2+1個點)(時域采集N個點,頻域只需要觀察N/2+1個點)
2、MATLAB中FFT的頻譜,應該看幅值
3、X軸頻率點的設置:采樣頻率為Fs,頻譜圖顯示的最高頻率為Fs/2(采樣定理)
:X軸頻率點:(0:1:N/2)*Fs/N
4、復數幅值修正
5、
/****************************************************/
/****************************************************/
/****************************************************/
栗子及實踐部分
一、信號
%% FFT clear;clc;close all Fs=1000; % 采集頻率 T=1/Fs; % 采集時間間隔 N=2000; % 采集信號的長度--采樣點數 f1=33; % 第一個余弦信號的頻率 f2=200; % 第二個余弦信號的頻率 t=(0:1:N-1)*T; % 定義整個采集時間點 t=t'; % 轉置成列向量 y=1.2+2.7*cos(2*pi*f1*t+pi/4)+5*cos(2*pi*f2*t+pi/6); % 時域信號
二、繪制時域信號
%% 繪制時域信號 figure plot(t,y) xlabel('時間') ylabel('信號值') title('時域信號')
三、FFT變換、並繪制-幅值、實部、虛部
%% fft變換 Y=fft(y); % Y為fft變換的結果,為復數向量 A=abs(Y); % 復數的幅值(模) RE=real(Y); % 復數的實部 IM=imag(Y); % 復數的虛部 %% 繪制fft變換結果(幅值,實部,虛部) figure subplot(3,1,1) plot(0:1:N-1,A) xlabel('序號 0 ~ N-1') ylabel('幅值') grid on %% 頻域只讀取0:1:N/2 subplot(3,1,2) plot(0:1:N-1,RE) xlabel('序號 0 ~ N-1') ylabel('實部') grid on subplot(3,1,3) plot(0:1:N-1,IM) xlabel('序號 0 ~ N-1') ylabel('虛部') grid on
可以看出頻域中的點關於(N/2)對稱,所以只需要觀察(0:1:N/2)
四、更改相位
幅值不受影響,但實部或虛部的值,會出現0的情況==>看MATLAB中FFT的頻譜,應該看幅值
繪制半譜圖(幅值的)后--我們發現-幅值-相位-頻率---均和時域對應不上。
==>進行幅值-修正
五、進行幅值-修正--並繪制圖形
%% fft變換 Y=fft(y); % Y為fft變換結果,復數向量 Y=Y(1:N/2+1); % 只看變換結果的一半即可 A=abs(Y); % 復數的幅值(模) f=(0:1:N/2)*Fs/N; % 生成頻率范圍 f=f'; % 轉置成列向量 %% 幅值修正 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函數的返回結果為弧度 phase_angle=rad2deg(phase_angle); plot(f,phase_angle) xlabel('頻率 (Hz)') ylabel('相位角 (degree)') title('FFT變換相位圖') grid on
放大后可以看到,此時,幅值-頻率都和時域一致
此時FFT的相位圖是雜亂無章的--不用擔心,沒有頻率處的相位是無意義的--我們只需要放大看各個(實際存在的)頻率點的相位即可
可以看到--f1=33Hz處為45度,即pi/4--是正確的
六、實際操作:請分析一個未知的采集信號 (example.mat),並確定該采集信號的頻率成分。其中, 信號的采集頻率 Fs = 2500 Hz
代碼
clear;clc;close all load('example') Fs=2500; % 采集頻率 T=1/Fs; % 采集時間間隔 N=length(y); % 采集信號的長度 t=(0:1:N-1)*T; % 定義整個采集時間點 t=t'; % 轉置成列向量 % 繪制時域信號 figure plot(t,y) xlabel('時間') ylabel('信號值') title('時域信號') % fft變換 Y=fft(y); % Y為fft變換結果,復數向量 Y=Y(1:N/2+1); % 只看變換結果的一半即可 A=abs(Y); % 復數的幅值(模) f=(0:1:N/2)*Fs/N; % 生成頻率范圍 f=f'; % 轉置成列向量 % 幅值修正 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函數的返回結果為弧度 phase_angle=rad2deg(phase_angle); plot(f,phase_angle) xlabel('頻率 (Hz)') ylabel('相位角 (degree)') title('FFT變換相位圖') grid on