1, Get code method

Get code method 1:
The complete code has been uploaded to my resources: [digital signal processing] synchronous compression and transformation of digital signals based on MATLAB [including Matlab source code 1535]

Get code method 2:
By subscribing to the payment column of zijishenguang blog, private bloggers can obtain this code with payment vouchers.

Note: if you subscribe to the paid column of zijishenguang blog, you can get a code for free (valid for three days from the Subscription Date);

2, Introduction to synchronous compression transformation

Based on WT, SST uses synchronous compression operator to improve the resolution of time-frequency ridge in time spectrum, and realizes the extraction and reconstruction of instantaneous frequency. set up ψ (b) Is the wavelet generating function, then the continuous wavelet transform of signal x(t) is:

Where:
x(t) - vibration signal;
W(a,b) - continuous wavelet transform result of x(t);
t -- time variable;
a - scale factor;
b -- translation factor;

According to the analysis, the instantaneous frequency information of position (a,b) in the wavelet domain is:

Where:
ω x(a,b) - instantaneous frequency;
j -- imaginary unit.
Literature [10] found that no matter what value a is, the oscillation characteristics of W(a,b) on B point to the initial frequency Ω. Therefore:

According to the defined synchronous compression transform, the inverse wavelet transform is:

Where:
x(b) - result of inverse wavelet transform;
C ψ—— Phase difference coefficient;
ψ (a ξ)—— Wavelet generating function.
yes ω x(a,b) is integrated along the direction of scale a and classified into frequency domain ω=ω At the position of x(a,b), the synchronous compression transformation is defined as:

Where:
Sst ( ω, b) - synchronous compression function of signal B;
ω—— Angular frequency.
The result of equation (6) and phase difference coefficient C ψ, The amplitude of the signal is reduced to the position in the frequency domain, and finally the high-resolution time spectrum is obtained.

3, Partial source code

% 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;

4, Operation results

5, matlab version and references

1 matlab version
2014a

2 references
[1] Shen Zaiyang. Proficient in MATLAB signal processing [M]. Tsinghua University Press, 2015
[2] Gao Baojian, Peng Jinye, Wang Lin, pan Jianshou. Signal and system -- Analysis and implementation using MATLAB [M]. Tsinghua University Press, 2020
[3] Wang Wenguang, Wei Shaoming, Ren Xin. MATLAB implementation of signal processing and system analysis [M]. Electronic Industry Press, 2018
[4] Gao Yanyan, Zhang Jing, Li Li, Jia Yingxi. Design of teaching demonstration system of digital signal processing based on GUI [J]. Education and Teaching Forum. 2019, (48)
[5] Li Jun, Zhang Shuling, Shuai Jing. Digital signal processing aided teaching system based on Matlab GUI interface [J]. Information communication. 2020, (08)
[6] Wang Bing, Wei Zhiheng, Wang Wenbin, Dai Yuanting, Zhao Minamata Jun. application of order analysis method based on synchronous compression transformation in bearing fault diagnosis of urban rail transit vehicles [J]. Research on urban rail transit. 2021,24 (07)