【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1534期】

原创

佛怒火莲 2021-12-01 14:37:21 博主文章分类：matlab完整代码 ©著作权

文章标签 matlab 自动驾驶开发语言数字信号处理参考文献 文章分类 代码人生

©著作权归作者所有：来自51CTO博客作者佛怒火莲的原创作品，请联系作者获取转载授权，否则将追究法律责任

二、同步压缩变换简介

SST以WT为基础，利用同步压缩算子提高时频脊线在时频谱上的分辨率，实现了对瞬时频率的提取与重构。设ψ(b)为小波母函数，则信号x(t)的连续小波变换为：

【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1534期】_参考文献

式中：

x(t)——振动信号；

W(a,b)——x(t)的连续小波变换结果；

t——时间变量；

a——尺度因子；

b——平移因子；

【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1534期】_参考文献_02

通过分析可知，在小波域中(a,b)位置的瞬时频率信息为：

【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1534期】_参考文献_03

式中：

ωx(a,b)——瞬时频率；

j——虚数单位。

文献[10]发现，无论a取何值，W(a,b)在b上的振荡特性均指向初始频率Ω,因此：

【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1534期】_自动驾驶_04

根据定义的同步压缩变换，小波逆变换为：

【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1534期】_开发语言_05

式中：

x(b)——小波逆变换结果；

Cψ——相差系数；

ψ(aξ)——小波母函数。

对ωx(a,b)沿尺度a方向整合，归划到频域中ω=ωx(a,b)的位置上，则同步压缩变换定义为：

【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1534期】_自动驾驶_06

式中：

Sst (ω,b)——信号b的同步压缩函数；

ω——角频率。

通过式(6)的结果与相差系数Cψ,将信号的幅值归化到时频域中所在位置，最终获得高分辨率的时频谱。

三、部分源代码

% A numerical signal.
clear;
SampFreq = 100;

t = 0 : 1/SampFreq : 14-1/SampFreq;


Sig = [sin(2*pi*(25*t))];

[m,n]=size(Sig);

time=(1:n)/SampFreq;
fre=(SampFreq/2)/(n/2):(SampFreq/2)/(n/2):(SampFreq/2);


Ts  = SST(Sig',100);

figure
imagesc(time,fre,abs(Ts));
axis xy
ylabel('Freq / Hz');
xlabel('Time / Sec');
title('SST');

%signal reconstruction.
s=real(sum(Ts));

%Reconstructed signal.
figure
plot(s);
title('Reconstructed signal');

%Original signal
figure
plot(Sig,'r-');
title('Original signal');
function [Ts] = SST(x,hlength);
% Computes the SST (Ts)  of the signal x.
% INPUT
%    x      :  Signal needed to be column vector.
%    hlength:  The hlength of window function.
% OUTPUT
%    Ts     :  The SST

[xrow,xcol] = size(x);

if (xcol~=1),
 error('X must be column vector');
end; 

if (nargin < 1),
error('At least 1 parameter is required');
end;

if (nargin < 2),
hlength=round(xrow/5);
end;

Siglength=xrow;
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';

% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
% derivative of window
dh = -2*pi/0.32^2*ht .* h; % g'

[hrow,hcol]=size(h); Lh=(hrow-1)/2; 

N=xrow;
t=1:xrow;

[trow,tcol] = size(t);


tfr1= zeros (N,tcol) ; 
tfr2= zeros (N,tcol) ; 

tfr= zeros (round(N/2),tcol) ; 
Ts= zeros (round(N/2),tcol) ;