如何使用matlab進行頻域分析

來自專欄 MATLAB學習12 人贊了文章

Matlab可以說是一個非常有用且功能齊全的工具，在通信、自控、金融等方面有廣泛的應用。

本文討論使用Matlab對信號進行頻域分析的方法。

說到頻域，不可避免的會提到傅里葉變換，傅里葉變換提供了一個將信號從時域轉變到頻域的方法。之所以要有信號的頻域分析，是因為很多信號在時域不明顯的特徵可以在頻域下得到很好的展現，可以更加容易的進行分析和處理。

FFT

Matlab提供的傅里葉變換的函數是FFT，中文名叫做快速傅里葉變換。快速傅里葉變換的提出是偉大的，使得處理器處理數字信號的能力大大提升，也使我們生活向數字化邁了一大步。

接下來就談談如何使用這個函數。

fft使用很簡單，但是一般信號都有x和y兩個向量，而fft只會處理y向量，所以想讓頻域分析變得有意義，那麼就需要用戶自己處理x向量

一個簡單的例子

從一個簡單正弦信號開始吧，正弦信號定義為：

我們現在通過以下代碼在Matlab中畫出這個正弦曲線

fo = 4; %frequency of the sine waveFs = 100; %sampling rateTs = 1/Fs; %sampling time intervalt = 0:Ts:1-Ts; %sampling periodn = length(t); %number of samplesy = 2*sin(2*pi*fo*t); %the sine curve %plot the cosine curve in the time domainsinePlot = figure;plot(t,y)xlabel(time (seconds))ylabel(y(t))title(Sample Sine Wave)grid

這就是我們得到的：

當我們對這條曲線fft時，我們希望在頻域得到以下頻譜（基於傅里葉變換理論，我們希望看見一個幅值為1的峰值在-4Hz處，另一個在+4Hz處）

使用FFT命令

我們知道目標是什麼了，那麼現在使用Matlab的內建的FFT函數來重新生成頻譜

%plot the frequency spectrum using the MATLAB fft commandmatlabFFT = figure; %create a new figureYfreqDomain = fft(y); %take the fft of our sin wave, y(t)stem(abs(YfreqDomain)); %use abs command to get the magnitude%similary, we would use angle command to get the phase plot!%well discuss phase in another post though!xlabel(Sample Number)ylabel(Amplitude)title(Using the Matlab fft command)gridaxis([0,100,0,120])

但是注意一下，這並不是我們真正想要的，有一些信息是缺失的

• x軸本來應該給我們提供頻率信息，但是你能讀出頻率嗎？
• 幅度都是100
• 沒有讓頻譜中心為

為FFT定義一個函數來獲取雙邊頻譜

以下代碼可以簡化獲取雙邊頻譜的過程，複製並保存到你的.m文件中

function [X,freq]=centeredFFT(x,Fs)%this is a custom function that helps in plotting the two-sided spectrum%x is the signal that is to be transformed%Fs is the sampling rateN=length(x);%this part of the code generates that frequency axisif mod(N,2)==0k=-N/2:N/2-1; % N evenelsek=-(N-1)/2:(N-1)/2; % N oddendT=N/Fs;freq=k/T; %the frequency axis%takes the fft of the signal, and adjusts the amplitude accordinglyX=fft(x)/N; % normalize the dataX=fftshift(X); %shifts the fft data so that it is centered

這個函數輸出正確的頻域範圍和變換後的信號，它需要輸入需要變換的信號和採樣率。

接下來使用前文的正弦信號做一個簡單的示例，注意你的示例.m文件要和centeredFFT.m文件在一個目錄下

[YfreqDomain,frequencyRange] = centeredFFT(y,Fs);centeredFFT = figure;%remember to take the abs of YfreqDomain to get the magnitude!stem(frequencyRange,abs(YfreqDomain));xlabel(Freq (Hz))ylabel(Amplitude)title(Using the centeredFFT function)gridaxis([-6,6,0,1.5])

效果如下：

這張圖就滿足了我們的需求,我們得到了在+4和-4處的峰值，而且幅值為1.

為FFT定義一個函數來獲取右邊頻譜

從上圖可以看出，FFT變換得到的頻譜是左右對稱的，因此，我們只需要其中一邊就能獲得信號的所有信息，我們一般保留正頻率一側。

以下的函數對上面的自定義函數做了一些修改，讓它可以幫助我們只畫出信號的正頻率一側

function [X,freq]=positiveFFT(x,Fs)N=length(x); %get the number of pointsk=0:N-1; %create a vector from 0 to N-1T=N/Fs; %get the frequency intervalfreq=k/T; %create the frequency rangeX=fft(x)/N; % normalize the data%only want the first half of the FFT, since it is redundantcutOff = ceil(N/2);%take only the first half of the spectrumX = X(1:cutOff);freq = freq(1:cutOff);

和前面一樣，使用正弦信號做一個示例，下面是示例代碼

[YfreqDomain,frequencyRange] = positiveFFT(y,Fs);positiveFFT = figure;stem(frequencyRange,abs(YfreqDomain));set(positiveFFT,Position,[500,500,500,300])xlabel(Freq (Hz))ylabel(Amplitude)title(Using the positiveFFT function)gridaxis([0,20,0,1.5])

效果如下：

--------------------------------------------

本文內容主要來自Matlab官方FFT教程，但是我找不到原鏈接了，就把我電腦上的文件上傳到了百度雲盤，公眾號後台回復

下載|FFT_tutorial

獲取下載鏈接

------------------------------------------------

本文代碼以上傳github，可自行下載測試，鏈接：

https://github.com/greedyhao/DSP/tree/master/matlab_tutorial

-------------------------------------------------