使用Matlab對(duì)信號(hào)進(jìn)行頻域分析的方法
Matlab可以說(shuō)是一個(gè)非常有用且功能齊全的工具,在通信、自控、金融等方面有廣泛的應(yīng)用。
本文引用地址:http://m.butianyuan.cn/article/201808/391336.htm本文討論使用Matlab對(duì)信號(hào)進(jìn)行頻域分析的方法。
說(shuō)到頻域,不可避免的會(huì)提到傅里葉變換,傅里葉變換提供了一個(gè)將信號(hào)從時(shí)域轉(zhuǎn)變到頻域的方法。之所以要有信號(hào)的頻域分析,是因?yàn)楹芏嘈盘?hào)在時(shí)域不明顯的特征可以在頻域下得到很好的展現(xiàn),可以更加容易的進(jìn)行分析和處理。
FFT
Matlab提供的傅里葉變換的函數(shù)是FFT,中文名叫做快速傅里葉變換??焖俑道锶~變換的提出是偉大的,使得處理器處理數(shù)字信號(hào)的能力大大提升,也使我們生活向數(shù)字化邁了一大步。
接下來(lái)就談?wù)勅绾问褂眠@個(gè)函數(shù)。
fft使用很簡(jiǎn)單,但是一般信號(hào)都有x和y兩個(gè)向量,而fft只會(huì)處理y向量,所以想讓頻域分析變得有意義,那么就需要用戶自己處理x向量
一個(gè)簡(jiǎn)單的例子
從一個(gè)簡(jiǎn)單正弦信號(hào)開(kāi)始吧,正弦信號(hào)定義為:
我們現(xiàn)在通過(guò)以下代碼在Matlab中畫(huà)出這個(gè)正弦曲線
fo = 4; %frequency of the sine wave
Fs = 100; %sampling rate
Ts = 1/Fs; %sampling time interval
t = 0:Ts:1-Ts; %sampling period
n = length(t); %number of samples
y = 2*sin(2*pi*fo*t); %the sine curve
%plot the cosine curve in the time domain
sinePlot = figure;
plot(t,y)
xlabel('time (seconds)')
ylabel('y(t)')
title('Sample Sine Wave')
grid
這就是我們得到的:
當(dāng)我們對(duì)這條曲線fft時(shí),我們希望在頻域得到以下頻譜(基于傅里葉變換理論,我們希望看見(jiàn)一個(gè)幅值為1的峰值在-4Hz處,另一個(gè)在+4Hz處)
使用FFT命令
我們知道目標(biāo)是什么了,那么現(xiàn)在使用Matlab的內(nèi)建的FFT函數(shù)來(lái)重新生成頻譜
%plot the frequency spectrum using the MATLAB fft command
matlabFFT = figure; %create a new figure
YfreqDomain = 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!
%we'll discuss phase in another post though!
xlabel('Sample Number')
ylabel('Amplitude')
title('Using the Matlab fft command')
grid
axis([0,100,0,120])
效果如下:
但是注意一下,這并不是我們真正想要的,有一些信息是缺失的
x軸本來(lái)應(yīng)該給我們提供頻率信息,但是你能讀出頻率嗎?
幅度都是100
沒(méi)有讓頻譜中心為
為FFT定義一個(gè)函數(shù)來(lái)獲取雙邊頻譜
以下代碼可以簡(jiǎn)化獲取雙邊頻譜的過(guò)程,復(fù)制并保存到你的.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 rate
N=length(x);
%this part of the code generates that frequency axis
if mod(N,2)==0
k=-N/2:N/2-1; % N even
else
k=-(N-1)/2:(N-1)/2; % N odd
end
T=N/Fs;
freq=k/T; %the frequency axis
%takes the fft of the signal, and adjusts the amplitude accordingly
X=fft(x)/N; % normalize the data
X=fftshift(X); %shifts the fft data so that it is centered
這個(gè)函數(shù)輸出正確的頻域范圍和變換后的信號(hào),它需要輸入需要變換的信號(hào)和采樣率。
接下來(lái)使用前文的正弦信號(hào)做一個(gè)簡(jiǎn)單的示例,注意你的示例.m文件要和centeredFFT.m文件在一個(gè)目錄下
[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')
grid
axis([-6,6,0,1.5])
效果如下:
這張圖就滿足了我們的需求,我們得到了在+4和-4處的峰值,而且幅值為1.
為FFT定義一個(gè)函數(shù)來(lái)獲取右邊頻譜
從上圖可以看出,F(xiàn)FT變換得到的頻譜是左右對(duì)稱(chēng)的,因此,我們只需要其中一邊就能獲得信號(hào)的所有信息,我們一般保留正頻率一側(cè)。
以下的函數(shù)對(duì)上面的自定義函數(shù)做了一些修改,讓它可以幫助我們只畫(huà)出信號(hào)的正頻率一側(cè)
function [X,freq]=positiveFFT(x,Fs)
N=length(x); %get the number of points
k=0:N-1; %create a vector from 0 to N-1
T=N/Fs; %get the frequency interval
freq=k/T; %create the frequency range
X=fft(x)/N; % normalize the data
%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2);
%take only the first half of the spectrum
X = X(1:cutOff);
freq = freq(1:cutOff);
和前面一樣,使用正弦信號(hào)做一個(gè)示例,下面是示例代碼
[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')
grid
axis([0,20,0,1.5])
效果如下:
評(píng)論