【信号处理】基于MATLAB短时傅里叶变换和小波变换的时频分析_聚沙成塔

[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。

MATLAB程序：

[plain]view plain copy

  [Au, Fs]=audioread('C:\Users\CDQ\Desktop\output.mp3');   % Fs 采样率 44100  

    [B, F, T, P] = spectrogram(Au(:,1),1024,512,1024,Fs);   % B是F大小行T大小列的频率峰值，P是对应的能量谱密度  

    figure  

    imagesc(T,F,C);  

    set(gca,'YDir','normal')  

    colorbar;  

    xlabel('时间 t/s');  

    ylabel('频率 f/Hz');  

    title('短时傅里叶时频图');

运行结果：
http://img.blog.csdn.net/20161204174629638?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center

2、小波变换

首先，在matlab中，小波变换的分析函数为cwt，其使用情况如下：

cwt函数功能：实现一维连续小波变换的函数。
cwt函数语法格式：

COEFS=cwt(S, SCALES, 'wname')

COEFS=cwt(S, SCALES, 'wname', 'plot')

COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE')

COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM)

使用说明：cwt为一维小波变换的函数。

格式 COEFS=cwt(S, SCALES, 'wname') 采用'wname'小波，在正、实尺度SCALES下计算向量一维小波系数。

格式 COEFS=cwt(S, SCALES, 'wname', 'plot') 除了计算小波系数外，还加以图形显示。

格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 计算并画出连续小波变换的系数，并使用PLOTMODE对图形着色。

格式 COEFS=cwt(S, SCALES, 'wname', 'plot') 相当于格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 中的语法 COEFS=cwt(S, SCALES, 'wname', 'absglb')

格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM) 能够计算并画出连续小波变换的系数。系数使用PLOTMODE和XLIM进行着色。其中：XLIM=[x1,x2]，并且有如下关系：1<=x1<=x2<=length(S)。

MODE值含义：

'lvl' scale-by-scale着色模式

'glb' 考虑所有尺度的着色模式

'abslvl'或'lvlabs' 使用系数绝对值的scale-by-scale着色模式

'absglb'或'glbabs' 使用系数绝对值并考虑所有尺度的着色模式

COEFS行的大小等于SCALES尺度的长度，COEFS列的大小等于信号S的长度。

MATLAB程序：

[plain]view plain copy

  wavename='cmor3-3';  

    totalscal=1024;  

    Fc=centfrq(wavename); % 小波的中心频率  测得Fc = 3  

    c=2*Fc*totalscal;    % 测得c = 1536  

    scals=c./(1:totalscal);  

    f=scal2frq(scals,wavename,1/fs); % 将尺度转换为频率   频率在0-500Hz取1024个点  

    coefs = cwt(Au(:,1),scals,wavename); % 求连续小波系数  

    t=0:1/fs:size(Au(:,1))/fs;  

    figure  

    imagesc(t,f,abs(coefs));  

    set(gca,'YDir','normal')  

    colorbar;  

    xlabel('时间 t/s');  

    ylabel('频率 f/Hz');  

    itle('小波时频图');

运行结果：

http://img.blog.csdn.net/20161204174640154?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：[转载]利用matlab怎样进行频谱分析

后一篇：【信号处理】滤波器的基本知识

新浪BLOG意见反馈留言板　欢迎批评指正