时频分析和MATLAB中的实现

来源：互联网发布：js 数组 clear 编辑：程序博客网时间：2024/06/06 03:32

文章来自

1. http://blog.csdn.net/hechongyang123/article/details/8501529

2. http://blog.sina.com.cn/s/blog_6163bdeb0102dwfw.html

在传统的信号处理中，人们分析和处理信号的最常用也是最直接的方法是傅里叶变换。傅里叶变换及其反变换构建起信号时域与频域之间变换的桥梁，是信号时域与频域分析的基础。但是以傅里叶变换为基础的经典分析方法，只是一种信号的整体变换，要么完全在时域进行，要么完全在频域进行，因而不具备时间和频率的定位功能，显然这对于平稳信号分析还是足够的。而对于非平稳信号而言，由于其频谱随时间有较大的变化，要求分析方法能够准确地反映出信号的局部时变频谱特性，只了解信号在时域或频域的全局特性是远远不够的，或者说是不适合的^[1]。例如，在地震勘探过程中，由于地层吸收等各种原因地震信号往往是非线性非平稳的。常规傅里叶变换和功率谱估计的方法就失去了物理意义，因此需要把整体谱推广到局部谱中来^[2]。时频分析方法是将一维时域信号映射到二维时频平面，全面地反映信号的时频联合特征。其基本思想是设计时间和频率的联合函数，以同时描述信号在不同时间和频率的能量密度和强度。如果有这样一个分布，就可以求在某一确定的频率和时间范围内的能量百分比，计算某一特定时间的频率密度及该分布的整体和局部的各阶矩等等。然而，不确定性原理不允许有“某个特定时间和频率处的能量”这一概念，因而理想的并不存在，只能研究伪能量密度或时频结构，根据不同的要求和不同的性能去逼近理想的时频表示^[3]。

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)|²，其中的的k是实值标量，定义如下
        对于单边PSD，计算公式如下，其中w(n)表示窗函数，Fs为采样频率，在0频率和奈奎斯特
        频率处，分子上的因子2改为1；

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

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