时频分析之短时傅里叶变换(STFT) |
您所在的位置：网站首页 › 傅里叶公式表 › 时频分析之短时傅里叶变换(STFT) |
时频分析之短时傅里叶变换(STFT)
|
目录 一、STFT 1.基本理论 2.spectrogram函数 3.频率分辨率和时间分辨率 3.1分辨率的影响因素 3.2提高频率分辨率的方法 二、MATLAB代码 参考文献 一、STFT 1.基本理论傅里叶变换只反映出信号在频域的特性，无法在时域内对信号进行分析。为了将时域和频域相联系，Gabor于1946年提出了短时傅里叶变换(short-time Fourier transform ，STFT)，其实质是加窗的傅里叶变换。STFT的过程是：在信号做傅里叶变换之前乘一个时间有限的窗函数 h(t)，并假定非平稳信号在分析窗的短时间隔内是平稳的，通过窗函数 h(t)在时间轴上的移动，对信号进行逐段分析得到信号的一组局部“频谱”。信号
由上式知，信号 短时傅里叶变换的优缺点：优点是其基本算法即是傅里叶变换，易于解释其物理意义；缺点是STFT 的窗宽是固定的，不能进行自适应调整。 2.spectrogram函数function varargout = spectrogram(x,varargin)：用于短时傅里叶变换（STFT）的频谱图。 S = SPECTROGRAM(X)，返回由矩阵S中的向量X指定的信号的短时傅立叶变换。默认情况下，X被分成8个有50个重叠的段，每个段用Hamming窗口加窗。用于计算离散傅里叶变换的频率点的数目等于256或大于段长度的下一次幂的256。如果X不能精确地分成八段，X将被截断。 S = SPECTROGRAM(X,WINDOW)，当WINDOW是一个向量时，将X分成与WINDOW长度相同的段，然后用WINDOW指定的向量对每个段进行加窗。如果WINDOW是一个整数，函数将X分成长度等于该整数值的段，并用Hamming窗口对每个段进行加窗。如果未指定WINDOW，则使用默认值。 S = SPECTROGRAM(X,WINDOW,NOVERLAP)，指定相邻线段之间重叠的NOVERLAP采样。如果WINDOW是整数，则NOVERLAP必须是小于WINDOW的整数。如果WINDOW是向量，则NOVERLAP必须是小于WINDOW长度的整数。如果未指定NOVERLAP，则使用默认值获得50个重叠。 S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT)，指定用于计算离散傅里叶变换的频率点的数目。如果未指定NFFT，则使用默认NFFT。 S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT,Fs)，指定采样率Fs（Hz）。如果Fs指定为空，则默认为1Hz。如果未指定，则使用规范化频率。 S的每一列包含X的短期、时间局部化频率含量的估计值。S列中的时间从左到右递增。频率从0开始沿行递增。如果X是一个长度为NX的复信号，则S是一个具有NFFT行和k=fix（（NX-NOVERLAP）/（length（WINDOW）-NOVERLAP））列的复矩阵。对于实X，如果NFFT是偶数，S有（NFFT/2+1）行；如果NFFT是奇数，S有（NFFT+1）/2行。 [S,F,T] = SPECTROGRAM(...)，返回频率向量F和时间向量T，在该向量下计算频谱图。F的长度等于S的行数。T的长度为k（如上定义），其值对应于每个线段的中心。如果没有提供采样率，F包含归一化频率。 [S,F,T] = SPECTROGRAM(X,WINDOW,NOVERLAP,F)，计算矢量F中指定的归一化频率下的双面谱图。F必须至少有两个元素。 [S,F,T] = SPECTROGRAM(X,WINDOW,NOVERLAP,F,Fs)，计算矢量F中指定频率下的双侧频谱图。F必须以赫兹表示，并且至少有两个元素。 [S,F,T,P] = SPECTROGRAM(...)，P是表示每段功率谱密度（PSD）的矩阵。对于实信号，SPECTROGRAM返回每段PSD的单边修正周期图估计值；对于复信号，如果指定了频率向量，则返回双边PSD。 3.频率分辨率和时间分辨率 3.1分辨率的影响因素窗口的长度会影响时间分辨率和频率分辨率。窗口越长，截取的信号越长，时间分辨率越低，频率分辨率越高；窗口越短，截取的信号越短，时间分辨率越高，频率分辨率越低。 从定量的角度来看，STFT的时间分辨率取决于平移步长 利用spectrogram函数自带两种提高频率分辨率的方法 方法1：将小于某个数值的数归零，提高频率分辨率。 调用格式：[S,F,T,P] = SPECTROGRAM(...,'MinThreshold',THRESH)；当10*log10(P)的对应元素小于THRESH(以分贝为单位指定阈值)时，将P的元素设置为零。 方法2：利用功率谱将能量集中的原理，提高频率分辨率。 调用格式：[S,F,T,P] = SPECTROGRAM(...,'reassigned')；将每个PSD估计值重新设定为其重心位置。 二、MATLAB代码 clc; clear; close all; tic; %% 产生信号 fs = 1000; % 采样频率1KHz t = 0:1/fs:1-1/fs; N=size(t,2); f1=10; f2=30; x=3*cos(2*pi*f1*t)+5*sin(2*pi*f2*t); figure(1); plot(t,x); %% 短时傅里叶变换 wlen=500;%设置窗口长度。窗口越长时间分辨率越差，频率分辨率越好。 hop=1;%每次平移的步长，最小为1。越小图像时间精度越好，但计算量大。 x=wkeep1(x,N+1*wlen);%中间截断 h=hamming(wlen);%设置海明窗的窗长 [B, F, T, P] = spectrogram(x,h,wlen-hop,N,fs); % B是F行T列的频率峰值，P是对应的能量谱密度 figure(2); imagesc(T,F,P); set(gca,'YDir','normal') colorbar; xlabel('时间 t/s'); ylabel('频率 f/Hz'); title('短时傅里叶时频图'); toc; 参考文献[1]宫宇新,何满潮,汪政红,等.岩石破坏声发射时频分析算法与瞬时频率前兆研究[J].岩石力学与工程学报,2013,32(4):787-799. [2]时频分析之STFT：短时傅里叶变换的原理与代码实现(非调用Matlab API). [3]matlab时频分析之短时傅里叶变换 spectrogram. |
【本文地址】