spectrogram函数做短时傅里叶分析

整理自: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,分辨率越高,相应的冗余就越多,计算量越大,但计算机只要能承受,问题不大。

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

时间: 2024-11-13 10:05:11

spectrogram函数做短时傅里叶分析的相关文章

Matlab_spectrogram_短时傅里叶分析_实现与讨论 [未完成]

在语音与音乐处理过程中,常用到短时傅里叶变换(Short Time Fourier Transformation, STFT).在一些学习路径中,STFT也是学习小波之前的预备知识.本文简单实现了 Matlab 中 Spectrogram 函数,在没有小波知识支撑下,讨论了参数的选择,以及分辨率相关的问题. 参考博客: 来源:    作者: 来源:    作者: 来源:    作者: 短时傅里叶变换简介 自己写一下对STFT的认识,不一定对. 短时傅里叶变换,其实还是傅里叶变换,只不过把一段长信

通过函数实现打印*号组成的直角三角形,函数要求传入行数即可。在main 方法中,通过用户输入得到行数,然后调用函数做打印。

#include <stdio.h> /* 1.通过函数实现打印*号组成的直角三角形,函数要求传入行数即可.在main方法中,通过用户输入得到行数,然后调用函数做打印.三角形样式:********************* */ int sanjiao(int hang){ int i; int j; for(i = 0; i < hang;i++) { for(j = 0;j <i+1;j++) { printf("*"); } printf("\n

JavaScript中使用函数做replace的第二个参数

通过第一个例子来全面看下作replace第二个参数的函数的参数....听起来是有点绕→_→ 例: //第一参数为正则表达式 1 var url = "http://www.softwhy.com/forum.php?mod=viewthread&tid=14743&extra=page%3D1"; 2 console.group("正则表达式"); 3 var regexp_global = /[?&](\w+)=([^&]*)/g;

UIApplicationMain函数做了什么?

在iOS应用中,每个程序得main函数中都调用了UIApplicationMain函数.   int main(int argc, char *argv[]) {     @autoreleasepool {        return UIApplicationMain(argc, argv, nil, NSStringFromClass([AppDelegate class]));     } } 先来看看UIApplicationMain函数的原型:   int UIApplication

RCurl getURL()函数做debug

getURL()函数做获取网页做debug,三步骤 1.首先创建一个对象debugGatherer(),该对象包含三个函数:(update(), value(), reset()); R> debugInfo <- debugGatherer() R> names(debugInfo) [1] "update" "value" "reset" R> class(debugInfo[[1]]) [1] "funct

用CIL写程序:写个函数做加法

前言: 上一篇文章小匹夫为CIL正名的篇幅比较多,反而忽略了写那篇文章初衷--即通过写CIL代码来熟悉它,了解它.那么既然有上一篇文章做基础(炮灰),想必各位对CIL的存在也就释然了,兴许也燃起了一点探索它,掌握它的欲望.那么小匹夫就继续扯一扯CIL,接下来的几篇文章也都以上一篇文章中的那个CIL实现的Hello Wolrd程序为基础,继续通过写CIL代码实现一些功能的方式来和各位探讨交流,同时也加深自己对CIL的掌握和印象. 人生就是做加法 "我的肩上搭着她得衣裳,我嗅着她留在衣服上的体香..

python中对函数做一些封装

把之前userlogin函数放到一个新的user.py中 在源文件中需要再做两步 1 import user 2 调用函数时userLogin 加上user. 变成user.userLogin

我们能用云函数做什么?

前言 本文以Firebase为例,因为腾讯云的云函数正在内测,还没申请到.:) 现如今云计算时代渐渐出现了越来越多的新型模式,从 IaaS: Infrastructure-as-a-Service(基础设施即服务) PaaS: Platform-as-a-Service(平台即服务) SaaS: Software-as-a-Service(软件即服务) 到CaaS:Containers as a Service(容器云) 再到的微服务架构,都在试着将各种软.硬件资源或抽象的事物做为一种服务提供给

2014.10.9 模板函数做友元函数的问题

1.模板函数作为类的友元函数,要先声名. 2.在类的定义中声明友元时要在函数后加上<T> 1 #ifndef Array1D_ 2 #define Array1D_ 3 4 // One-dimensional arrays. 5 6 #include <iostream> 7 #include "xcept.h" 8 9 10 template <class T> class Array1D; 11 //template <class T&g