码迷,mamicode.com
首页 > 其他好文 > 详细

spectrogram函数做短时傅里叶分析

时间:2019-01-02 12:44:59      阅读:188      评论:0      收藏:0      [点我收藏+]

标签:例子   频谱图   而在   oss   忽略   ...   surf   常量   list   

整理自:http://blog.sina.com.cn/s/blog_6163bdeb0102dwfw.html

今天偶人发现原来matlab自带了短时傅里叶变换的分析函数,老版本的matlab是specgram函数,新的改成了spectrogram函数,虽然一说到时频分析,都会说到小波分析,小波分析要比短时傅里叶要好云云,但在分析信号的瞬时频谱时,短时傅里叶还是有它的用武之地的。前一阵也看了一些有关小波分析的matlab实现,发现帮助中使用小波也多是除噪、压缩,都说小波是时频显微镜,它的用武之地还是在于查看高频在哪一级分解中,进而可以有效滤除一些信号,比如除噪,所以短时傅里叶变换查看瞬时频率正好互补一下。时频分析还认识的不深,一个阶段的想法而已。

另外,之前对matlab的扫频函数chirp做过总结,见http://blog.sina.com.cn/s/blog_6163bdeb0100qbqo.html,里面就是使用spectrogram函数来查看产生的扫频信号的瞬时频率的,当时不知道那个函数是干啥,就感觉好神奇,现在正好看到,总结一下吧!

 

spectrogram

功能:使用短时傅里叶变换得到信号的频谱图。

语法:

       [S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)

       [S,F,T,P]=spectrogram(x,window,noverlap,F,fs)

说明: 当使用时无输出参数,会自动绘制频谱图;有输出参数,则会返回输入信号的短时傅里叶变

         换。当然也可以从函数的返回值S,F,T,P绘制频谱图,具体参见例子。

参数:

x---输入信号的向量。默认情况下,即没有后续输入参数,x将被分成8段分别做变换处理,

        如果x不能被平分成8段,则会做截断处理。默认情况下,其他参数的默认值为

          window---窗函数,默认为nfft长度的海明窗Hamming

          noverlap---每一段的重叠样本数,默认值是在各段之间产生50%的重叠

          nfft---做FFT变换的长度,默认为256和大于每段长度的最小2次幂之间的最大值。

                   另外,此参数除了使用一个常量外,还可以指定一个频率向量F

          fs---采样频率,默认值归一化频率

Window---窗函数,如果window为一个整数,x将被分成window段,每段使用Hamming窗函数加窗。

         如果window是一个向量,x将被分成length(window)段,每一段使用window向量指定的

         窗函数加窗。所以如果想获取specgram函数的功能,只需指定一个256长度的Hann窗。

Noverlap---各段之间重叠的采样点数。它必须为一个小于window或length(window)的整数。

           其意思为两个相邻窗不是尾接着头的,而是两个窗有交集,有重叠的部分。

Nfft---计算离散傅里叶变换的点数。它需要为标量。

Fs---采样频率Hz,如果指定为[],默认为1Hz。

S---输入信号x的短时傅里叶变换。它的每一列包含一个短期局部时间的频率成分估计,

    时间沿列增加,频率沿行增加。

    如果x是长度为Nx的复信号,则S为nfft行k列的复矩阵,其中k取决于window,

        如果window为一个标量,则k = fix((Nx-noverlap)/(window-noverlap))

        如果window为向量,则k = fix((Nx-noverlap)/(length(window)-noverlap))

    对于实信号x,如果nfft为偶数,则S的行数为(nfft/2+1),如果nfft为奇数,

    则行数为(nfft+1)/2,列数同上。

F---在输入变量中使用F频率向量,函数会使用Goertzel方法计算在F指定的频率处计算频谱图。

    指定的频率被四舍五入到与信号分辨率相关的最近的DFT容器(bin)中。而在其他的使用nfft

    语法中,短时傅里叶变换方法将被使用。对于返回值中的F向量,为四舍五入的频率,其长度

    等于S的行数。

T---频谱图计算的时刻点,其长度等于上面定义的k,值为所分各段的中点。

P---能量谱密度PSD(Power Spectral Density)

    对于实信号,P是各段PSD的单边周期估计;

    对于复信号,当指定F频率向量时,P为双边PSD。

    P矩阵的元素计算公式如下P(I,j)=k|S(I,j)|2,其中的的k是实值标量,定义如下

        对于单边PSD,计算公式如下,其中w(n)表示窗函数,Fs为采样频率,在0频率和奈奎斯特

        频率处,分子上的因子2改为1;

                                技术分享图片

        对于双边PSD,计算公式如下

                               技术分享图片

         如果采样频率没有指定,分母上的Fs由2*pi代替。

 

spectrogram(...)当调用函数时没有输出参数,将会自动绘制各段的PSD估计,绘制的命令如下

       surf(T,F,10*log10(abs(P)));

       axis tight;

       view(0,90);

spectrogram(...,‘freqloc‘)使用freqloc字符串可以控制频率轴显示的位置。当freqloc=xaxis

       时,频率轴显示在x轴上,当freqloc=yaxis时,频率轴显示在y轴上,默认是显示在x轴

       上。如果在指定freqloc的同时,又有输出变量,则freqloc将被忽略。

 

例.计算并显示二次扫频信号的PSD图,扫频信号的频率开始于100Hz,在1s时经过200Hz

T = 0:0.001:2;

X = chirp(T,100,1,200,‘q‘);

spectrogram(X,128,120,128,1E3);

title(‘Quadratic Chirp‘);

技术分享图片

 频率显示在y轴上:

t=0:0.001:2; % 2 secs @ 1kHz sample rate
y=chirp(t,100,1,200,‘q‘); % Start @ 100Hz, cross 200Hz at t=1sec
spectrogram(y,kaiser(128,18),120,128,1E3,‘yaxis‘);
title(‘Quadratic Chirp: start at 100Hz and cross 200Hz at t=1sec‘);

技术分享图片

 

例.计算并显示线性扫频信号的PSD图,扫频信号由直流开始,在1s时经过150Hz,控制频率轴显示在y轴上

T = 0:0.001:2;

X = chirp(T,0,1,150);

[S,F,T,P] = spectrogram(X,256,250,256,1E3);

surf(T,F,10*log10(P),‘edgecolor‘,‘none‘); axis tight;

view(0,90);

xlabel(‘Time (Seconds)‘); ylabel(‘Hz‘);

技术分享图片

 

函数使用的注意:

nfft越大,频域的分辨率就越高(分辨率=fs/nfft),但离瞬时频率就越远;

noverlap影响时间轴的分辨率,越接近nfft,分辨率越高,相应的冗余就越多,计算量越大,但计算机只要能承受,问题不大。

spectrogram函数做短时傅里叶分析

标签:例子   频谱图   而在   oss   忽略   ...   surf   常量   list   

原文地址:https://www.cnblogs.com/shuqingstudy/p/10207354.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!