时频分析-短时傅里叶、小波变换(一)

目录：

1.三种时频方式解决的问题
2.短时傅里叶变换
3.小波变换

一、三种时频变换方式解决的问题

（一）傅里叶变换

FT在平稳信号的分析和处理中有着突出贡献的原因在于，人们利用它可以把复杂的时间信号和空间信号变换到频率域中，然后用相对简单的频谱特性去分析和发现原信号的动态特性。

FT 正变换告诉我们：从时间（空间）信号中提取信号的频谱信息F(W)，就是使用整个时间域的所有信息来计算单个确定频率的谱值（频域函数F(W)的任一频率w 对应的函数值），这是由时间轴（−∞,∞）上的确定信号f(t)决定的。因此，它求出的频域函数对应的时整个时间轴，所以可以知道，傅里叶变换对频谱的描绘是“全局性”的，不能反映时间维度局部区域上的特征，人们虽然从傅立叶变换能清楚地看到一整段信号包含的每一个频率的分量值，但很难看出对应于频率域成分的不同时间信号的持续时间和发射的持续时间，缺少时间信息使得傅立叶分析再更精密的分析中失去作用。

伊利诺依斯大学教授曾说：“若你记录1小时长的信息而在最后5分钟出错，这一错误就会毁了整个傅立叶变换。相位的错误是灾难性的，如果在相位上哪怕犯了一个错误，你最后就会发现你所干的事与最初的信号无关了。”

(二)短时傅里叶变换（窗式傅里叶变换）

1. 基本思想
STFT（short-time Fourier transform，短时傅里叶变换）是和傅里叶变换相关的一种数学变换，用以确定时变信号其局部区域正弦波的频率与相位。用途：与小波变换相似，经STFT处理后的信号具有时域和频域的局部化特性（见下图），可以借助其分析信号的视频特性。

局部平稳化-把长的非平稳随机过程看成是一系列短时随机平稳信号的叠加，短时性可通过在时间上加窗口函数实现（即截取一部分源数据）。通过该方法，人们至少可以说，无论发现了什么频率成分，它一定是发生在信号被截取的某个特定时间段内。
短时傅里叶变换是最常用的一种时频分析方法，它通过时间窗内的一段信号来表示某一时刻的信号特征。在短时傅里叶变换过程中，窗的长度决定频谱图的时间分辨率和频率分辨率，窗长越长，截取的信号越长，信号越长，傅里叶变换后频率分辨率越高，时间分辨率越差；相反，窗长越短，截取的信号就越短，频率分辨率越差，时间分辨率越好，也就是说短时傅里叶变换中，时间分辨率和频率分辨率之间两者不可兼得，必须根据具体需求进行取舍。
简单来说，短时傅里叶变换就是先把一个函数和窗函数进行相乘，然后再进行一维的傅里叶变换。并通过窗函数的滑动得到一系列的频谱函数，将这些结果依次开便得到一个二维的时频图。

2.数学公式

3.matlab函数
首先，在matlab中，短时傅里叶变换的分析函数为spectrogram，其使用情况如下：

 [S,F,T,P]=spectrogram(x,window,noverlap,nfft,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。

%% 变频信号的fft观察
close all
clear,clc
%% 定义一个变频信号的参数
PIx2 = 2 * pi;
Fs = 1500;
T = 1/Fs;
LengthOfSignal = 3000;
NFFT = 2^nextpow2(LengthOfSignal); %要计算的点数
% NFFT = 128;
t = (0:LengthOfSignal-1)*T;
amp1 = 2;
amp2 = 2;
amp3 = 1.5;
offset = 2;
freq1 = 100;
freq2 = 150;
freq3 = 300;
signal = zeros(1,length(t));
 
%% 定义信号
for temp = 1:LengthOfSignal
    if(temp <= LengthOfSignal/4)
        signal(temp) =offset + amp1 * sin(PIx2 * freq1 * t(temp));
    elseif(temp <= LengthOfSignal/2)
        signal(temp) =offset -1*amp2 * sin(PIx2 * freq2 * t(temp));
    elseif(temp <= 3*LengthOfSignal/4)
        signal(temp) =offset -1*amp3 * sin(PIx2 * freq3 * t(temp));
    else
        signal(temp) =offset + amp1 * sin(PIx2 * freq1 * t(temp));
    end
end
 
%% 绘制图形
subplot(311);
plot(t, signal);
grid on
title('signal of different frequency');
xlabel('time');
ylabel('amp');
 
%% fft运算
fMax = NFFT/2 + 1;
signalFFT = abs(fft(signal,NFFT));
% signalFFTShow = 2 * abs(fft(signal(1:fMax),NFFT)/LengthOfSignal);
signalFFTShow = 2 * signalFFT / LengthOfSignal;
signalFFTShow(1) = signalFFTShow(1)/2;
f = Fs/2*linspace(0,1,fMax);
subplot(312);
plot(f,signalFFTShow(1:fMax));
grid on
title('fft signal');
xlabel('frequency');
ylabel('amp');
 
%% surf测试三维图(不合理x(j),y(i),z(i,j)对应)
 
subplot(313)
spectrogram(signal,128,120,128,1.5e3,'yaxis'); %时频域图

(三)小波变换

1. 思想方法

由短时傅立叶变换对函数（信号）进行的分析，相当于用一个形状、大小和放大倍数相同的“放大镜”在时－频域平面上移动去观察某固定长度时间内的频率特性。这里的问题是：尽管窗式傅立叶变换能解决变换函数的时域局域化问题，但是，其窗口的大小和形状是固定的，即窗口没有自适应性。和短时傅里叶变换相比，小波变换有着窗口自适应的特点，即高频信号分辨率高（但是频率分辨率差），低频信号频率分辨率高（但是时间分辨率差），而在工程中常常比较关心低频的频率，关心高频出现的时间，所以近些年用途比较广泛。

2. 数学公式

3. matlab函数

  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的长度。

clc
clear all
close all
% 原始信号
fs=1000;
f1=50;
f2=100;
t=0:1/fs:1;
s=sin(2*pi*f1*t)+sin(2*pi*f2*t);
figure
plot(t, s)
% 连续小波变换
wavename='cmor3-3';
totalscal=256;
Fc=centfrq(wavename); % 小波的中心频率
c=2*Fc*totalscal;
scals=c./(1:totalscal);
f=scal2frq(scals,wavename,1/fs); % 将尺度转换为频率
coefs=cwt(s,scals,wavename); % 求连续小波系数
figure
imagesc(t,f,abs(coefs));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('小波时频图');