原创 如何使用matlab进行频域分析

2018-8-24 20:08 10610 39 12 分类: 处理器与DSP

Matlab可以说是一个非常有用且功能齐全的工具,在通信、自控、金融等方面有广泛的应用。

本文讨论使用Matlab对信号进行频域分析的方法。

说到频域,不可避免的会提到傅里叶变换,傅里叶变换提供了一个将信号从时域转变到频域的方法。之所以要有信号的频域分析,是因为很多信号在时域不明显的特征可以在频域下得到很好的展现,可以更加容易的进行分析和处理。


FFT

Matlab提供的傅里叶变换的函数是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])


效果如下:



但是注意一下,这并不是我们真正想要的,有一些信息是缺失的

  • 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 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定义一个函数来获取右边频谱

    从上图可以看出,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

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

    个人博客


    文章评论2条评论)

    登录后参与讨论

    用户1208903 2018-9-5 12:19

    学习学习。

    Orima 2018-8-25 22:11

    好好参考...
    相关推荐阅读
    greedyhao 2019-03-17 09:03
    uboot链接脚本分析
    上一篇uboot主编译脚本分析的编译选项参数中提到了链接脚本这么个东西我也是写到这里才开始了解链接脚本(Linker command scripts),之前只是大概知道有这么一个东西先从参考资料中的内...
    greedyhao 2018-11-28 08:49
    51总线方式获取adc0809数值
    最近单片机有一个实验挺有意思的,使用51单片机以总线的方式读取adc0809的数据先补充点关于单片机总线的预备知识,我一开始不了解总线的时候做这个实验也是很懵逼的。单片机的三总线结构​51单片机有三条...
    greedyhao 2018-09-29 21:18
    对Xil_Out32未定义的引用
    第一次在HLS中遇见一个官方库函数未定义问题,这就把解决方法记录下来。在创建好工程,写完工程代码后,发现报错如下./src/led_controller_test_tut_4A.o:在函数‘main’...
    greedyhao 2018-09-29 21:12
    Vivado生成Bitstream失败的解决方法
    跟着实验指导书,难得的又遇到问题了,在最后生成Bitstream的时候出错了,无法生成Bitstream。报错信息如下[DRC NSTD-1] Unspecified I/O Standard: 4 ...
    greedyhao 2018-09-17 22:26
    zybo初体验
    前段时间过生日,亲戚给了些钱让我自己买生日礼物,早就想搞一块zynq的板了,但是淘宝一看太贵,于是去咸鱼上淘了一块二手的ZYBO(谁知道是不是二手,反正看着挺新)来玩玩。学习肯定是从官方资料开始的,T...
    greedyhao 2018-09-13 10:53
    运行c代码前发生了什么
    众所周知,c代码也是一种比较高级的语言了,机器是没有办法直接运行的,机器所能理解的只有机器码--那一串0和1而已。虽然早就知道c代码是先经过编译链接最后才放在机器上执行,但是在这么一个过程中究竟发生了...
    我要评论
    2
    39
    关闭 站长推荐上一条 /2 下一条