Matlab可以说是一个非常有用且功能齐全的工具,在通信、自控、金融等方面有广泛的应用。 本文讨论使用Matlab对信号进行频域分析的方法。 说到频域,不可避免的会提到傅里叶变换,傅里叶变换提供了一个将信号从时域转变到频域的方法。之所以要有信号的频域分析,是因为很多信号在时域不明显的特征可以在频域下得到很好的展现,可以更加容易的进行分析和处理。 FFTMatlab提供的傅里叶变换的函数是FFT,中文名叫做快速傅里叶变换。快速傅里叶变换的提出是伟大的,使得处理器处理数字信号的能力大大提升,也使我们生活向数字化迈了一大步。 接下来就谈谈如何使用这个函数。 fft使用很简单,但是一般信号都有x和y两个向量,而fft只会处理y向量,所以想让频域分析变得有意义,那么就需要用户自己处理x向量 一个简单的例子从一个简单正弦信号开始吧,正弦信号定义为: 我们现在通过以下代码在Matlab中画出这个正弦曲线 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 这就是我们得到的: 当我们对这条曲线fft时,我们希望在频域得到以下频谱(基于傅里叶变换理论,我们希望看见一个幅值为1的峰值在-4Hz处,另一个在+4Hz处) 使用FFT命令我们知道目标是什么了,那么现在使用Matlab的内建的FFT函数来重新生成频谱 %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])
但是注意一下,这并不是我们真正想要的,有一些信息是缺失的
以下代码可以简化获取双边频谱的过程,复制并保存到你的.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)/2N-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 这个函数输出正确的频域范围和变换后的信号,它需要输入需要变换的信号和采样率。 接下来使用前文的正弦信号做一个简单的示例,注意你的示例.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') grid axis([-6,6,0,1.5]) 效果如下: 这张图就满足了我们的需求,我们得到了在+4和-4处的峰值,而且幅值为1. 从上图可以看出,FFT变换得到的频谱是左右对称的,因此,我们只需要其中一边就能获得信号的所有信息,我们一般保留正频率一侧。 以下的函数对上面的自定义函数做了一些修改,让它可以帮助我们只画出信号的正频率一侧 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); 和前面一样,使用正弦信号做一个示例,下面是示例代码 [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]) 效果如下: ------------------------------------------------ 大部分内容是根据Matlab的FFT教程翻译的,但是源文档的下载链接找不到了 ------------------------------------------------ 本文代码以上传github,可自行下载测试,链接: https://github.com/greedyhao/DSP/tree/master/matlab_tutorial |
|