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

Posted 2021-12-22 紫极神光

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【数字信号处理】基于matlab数字信号同步压缩变换【含Matlab源码 1535期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、同步压缩变换简介

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

式中：
x(t)——振动信号；
W(a,b)——x(t)的连续小波变换结果；
t——时间变量；
a——尺度因子；
b——平移因子；

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

式中：
ωx(a,b)——瞬时频率；
j——虚数单位。
文献[10]发现，无论a取何值，W(a,b)在b上的振荡特性均指向初始频率Ω,因此：

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

式中：
x(b)——小波逆变换结果；
Cψ——相差系数；
ψ(aξ)——小波母函数。
对ωx(a,b)沿尺度a方向整合，归划到频域中ω=ωx(a,b)的位置上，则同步压缩变换定义为：

式中：
Sst (ω,b)——信号b的同步压缩函数；
ω——角频率。
通过式(6)的结果与相差系数Cψ,将信号的幅值归化到时频域中所在位置，最终获得高分辨率的时频谱。

三、部分源代码

% A numerical signal.
clear;
SampFreq = 100;
t1 = 0 : 1/SampFreq : 6;
t2 = 6+1/SampFreq : 1/SampFreq : 14-1/SampFreq;
t = [t1 t2];

Sig1 = [sin(2*pi*((25+8)*t1 + 10*sin(t1))) sin(2*pi*(34.2+8)*t2) ];
Sig2 = sin(2*pi*(8*t+3*atan((t - 5).^2)));
Sig=Sig1+Sig2;
[m,n]=size(Sig);
time=(1:n)/SampFreq;
fre=(SampFreq/2)/(n/2):(SampFreq/2)/(n/2):(SampFreq/2);

Ts  = SST(Sig',50);

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

df=fre(2)-fre(1);
%Signal Reconstrucion

%reconstruted region of Sig1 is 20-45Hz.
s1=real(sum(Ts(20/df:45/df,:)));
figure
plot(s1);hold on;
plot(Sig1,'r-');
title('Reconstructed signal (blue),original signal (red)');

%reconstruted region of Sig1 is 1-15Hz.
s2=real(sum(Ts(1/df:15/df,:)));
figure
plot(s2);hold on;
plot(Sig2,'r-');
title('Reconstructed signal (blue),original signal (red)');

%It provides a perfect reconstrution performance.
%Maybe you cannot see the bule signal clearly.
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) ; 

for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1; 
rSig = x(ti+tau,1);
%rSig = hilbert(real(rSig));
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr2(indices,icol)=rSig.*conj(dh(Lh+1+tau));
end;

tfr1=fft(tfr1);
tfr2=fft(tfr2);

tfr1=tfr1(1:round(N/2),:);
tfr2=tfr2(1:round(N/2),:);

ft = 1:round(N/2);
bt = 1:N;

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
[4]高彦彦,张晶,李莉,贾英茜.基于GUI的《数字信号处理》教学演示系统的设计[J].教育教学论坛. 2019,(48)
[5]李俊,张淑玲,帅晶.基于Matlab GUI界面的数字信号处理辅助教学系统[J].信息通信. 2020,(08)
[6]王冰,魏志恒,王文斌,戴源廷,赵俣钧.基于同步压缩变换的阶比分析法在城市轨道交通车辆轴承故障诊断中的应用[J].城市轨道交通研究. 2021,24(07)