FFT-Matlab初步實現


 

 

 

 

 

 

 

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

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

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

下面是具體說明

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

  

 

 

 

 

全部代碼下載地址

 

 

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM