0
点赞
收藏
分享

微信扫一扫

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

萨摩斯加士奇 2022-02-28 阅读 178

目录:

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

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

(一)傅里叶变换

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

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

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

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

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

2.数学公式

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

 [S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)
%% 变频信号的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函数

使用说明: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('小波时频图');

举报

相关推荐

0 条评论