同步壓縮變換原理
作為處理非平穩(wěn)信號(hào)的有力工具，時(shí)頻分析在時(shí)域和頻域聯(lián)合表征信號(hào)，是時(shí)間和頻率的二元函數(shù)。傳統(tǒng)的時(shí)頻分析工具主要分為線性方法和二次方法。線性方法受到海森堡測不準(zhǔn)原理的制約，二次方法存在交叉項(xiàng)的干擾。
為了提升時(shí)頻聚集性，逼近理想的時(shí)頻表示，時(shí)頻重排 (Reassignment method, RM)作為一種后處理技術(shù)被提。它在二維的時(shí)頻面上重排時(shí)頻系數(shù)，導(dǎo)致其喪失了重構(gòu)信號(hào)的能力。同步壓縮變換作為一種特殊的重排，不僅可以銳化時(shí)頻表示，還能恢復(fù)信號(hào)。因此，同步壓縮變換受到研究學(xué)者的熱愛。
同步壓縮變換主要有兩種框架，一種基于小波變換，另一種基于短時(shí)傅里葉變換。本文中以小波變換為框架，介紹同步壓縮變換。在2011年，Daubechies等人提出了同步壓縮小波變換，作為一種經(jīng)驗(yàn)?zāi)B(tài)分解工具。它的主要思想是將小波變換的系數(shù)重新分配到估計(jì)的瞬時(shí)頻率處。所以，同步壓縮變換的核心是估計(jì)瞬時(shí)頻率。通常，是對(duì)時(shí)頻表示關(guān)于時(shí)間求導(dǎo)，估計(jì)信號(hào)的瞬時(shí)頻率。然后，把系數(shù)擠壓到估計(jì)的瞬時(shí)頻率處。盡管，同步壓縮小波變換有很多改進(jìn)形式。但是，萬變不離其宗，都是同樣的套路。重排和同步壓縮變換的主要思想如下，
下圖是2011年論文中，一個(gè)諧波信號(hào)的小波變換和同步壓縮小波變換的結(jié)果對(duì)比。左邊是諧波信號(hào)的時(shí)域波形，中間是小波變換，右邊是同步壓縮小波變換的結(jié)果。顯然，小波變換的系數(shù)分布在瞬時(shí)頻率8Hz的周圍，同步壓縮小波變換的系數(shù)就只分布在8Hz處。換句話說，同步壓縮小波變換的時(shí)頻分辨率更高，能量聚集性更高。一般通過瑞利熵來衡量時(shí)頻表示的能量聚集性，但是對(duì)于真實(shí)信號(hào)這并不是合理的選擇。
同步壓縮小波變換代碼
主程序
clear ;
close all ;
clc;
%% fix the seed
initstate(23400) ;
opts = struct();
opts.motherwavelet = 'morlet' ;%morlet,Cinfc
opts.CENTER = 1;
opts.FWHM = 0.2;
fs=2^9; %采樣頻率
dt=1/fs; %時(shí)間精度
timestart=0;
timeend=1;
t=(0:(timeend-timestart)/dt-1)*dt+timestart;
t=t';
k1=35;
f1=60;
y=exp((1i)*pi*k1*t.*t+2*pi*1i*f1*t).*(t>=0& t<timeend);
freqlow=0;
freqhigh=128;
alpha=1;
seta=-acot(k1);
f=freqlow:alpha:freqhigh;
e=t(end);
%注意SQCWT中列為尺度，行為時(shí)間
[tfr, tfrsq, tfrtic, tfrsqtic] = sqCWT(t, y, freqlow, freqhigh, alpha, opts);
figure();
imagesc(t,tfrtic,abs(tfr));
set(gca,'ydir','normal');
xlabel('Time(s)','FontSize',12,'FontName','Times New Roman');
ylabel('Scales','FontSize',12,'FontName','Times New Roman');
set(gca,'FontSize',16,'FontName','Times New Roman')
set(gca,'linewidth',1);
%
figure()
imagesc(t,f, abs(tfrsq));
set(gca,'ydir','normal');
xlabel('Time(s)','FontSize',12,'FontName','Times New Roman');
ylabel('Frequency(Hz)','FontSize',12,'FontName','Times New Roman');
set(gca,'FontSize',16,'FontName','Times New Roman')
set(gca,'linewidth',1);
同步壓縮變換核心代碼
function [tfr, tfrsq, tfrtic, tfrsqtic] = sqCWT(t, x, freqlow, freqhigh, alpha, opts);
%
% Synchrosqueezing transform ver 0.5 (2015-03-09)
% You can find more information in
% http://sites.google.com/site/hautiengwu/
%
% Example: [~, tfrsq, ~, tfrsqtic] = sqCWT(time, xm, lowfreq, highfreq, alpha, opts);
% time: time of the signal
% xm: the signal to be analyzed
% [lowfreq, highfreq]: the frequency range in the output time-frequency representation. For the sake of computational efficiency.
% alpha: the frequency resolution in the output time-frequency representation
% opts: parameters for the CWT analysis. See below
% tfr/tfrtic: the CWT and its scale tic
% tfrsq/tfrsqtic: the SST-CWT and its frequency tic
%
% by Hau-tieng Wu v0.1 2011-06-20 (hauwu@math.princeton.edu)
% v0.2 2011-09-10
% v0.3 2012-03-03
% v0.4 2012-12-12
% v0.5 2015-03-09
%% you can play with these 4 parameters, but the results might not be
%% that different if the change is not crazy
Gamma = 1e-8;
nvoice = 32;
scale = 2;
oct = 1;
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
noctave = floor(log2(length(x))) - oct;
dt = t(2) - t(1);
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%% Continuous wavelet transform
[tfr, tfrtic] = CWT(t, x, oct, scale, nvoice, opts) ;
% tfr=tfr';
% tfrtic=tfrtic';
%行時(shí)間，列尺度
Dtfr = (-i/2/pi/dt)*[tfr(2:end,:) - tfr(1:end-1,:); tfr(end,:)-tfr(end-1,:)] ;
% Dtfr=(-1i/2/pi/dt)*[tfr(:,2:end)-tfr(:,1:end-1) tfr(:,end)-tfr(:,end-1)] /(dt);
Dtfr((abs(tfr) < Gamma)) = NaN;
omega = Dtfr./tfr;
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%% Synchro-squeezing transform
[tfrsq, tfrsqtic] = SQ(tfr, omega, alpha, scale, nvoice, freqlow, freqhigh);
tfr = transpose(tfr) ;
tfrsq = transpose(tfrsq) ;
end
%=====================================================================
%% function for CWT-based SST
function [tfrsq, freq] = SQ(tfd, omega, alpha, scale, nvoice, freqlow, freqhigh);
omega = abs(omega);
[n, nscale] = size(tfd);
nalpha = floor((freqhigh - freqlow)./alpha);
tfrsq = zeros(n, nalpha);
freq = ([1:1:nalpha])*alpha + freqlow;
nfreq = length(freq);
for b = 1:n %% Synchro-
for kscale = 1: nscale %% -Squeezing
% qscale = scale .* (2^(kscale/nvoice));
qscale = scale .* (2^(kscale/nvoice));
if (isfinite(omega(b, kscale)) && (omega(b, kscale)>0))
k = floor( ( omega(b,kscale) - freqlow )./ alpha )+1;
if (isfinite(k) && (k > 0) && (k < nfreq-1))
ha = freq(k+1)-freq(k);
tfrsq(b,k) = tfrsq(b,k) + log(2)*tfd(b,kscale)*sqrt(qscale)./ha/nvoice;
end
end
end
end
end
小波變換代碼
function [tfr, tfrtic] = newOmega(t, x, oct, scale, nvoice, opts)
%
% Continuous wavelet transform ver 0.1
% Modified from wavelab 85
%
% INPUT:
% OUTPUT:
% DEPENDENCY:
%
% by Hau-tieng Wu 2011-06-20 (hauwu@math.princeton.edu)
% by Hau-tieng Wu 2012-12-23 (hauwu@math.princeton.edu)
%
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%% prepare for the input
if nargin < 6
opts = struct();
opts.motherwavelet = 'Cinfc' ;
opts.CENTER = 1 ;
opts.FWHM = 0.3 ;
end
if nargin < 5
nvoice = 32;
end
if nargin < 4
scale = 2;
end
if nargin < 3
oct = 1;
end
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%% start to do CWT
dt = t(2)-t(1);
n = length(x);
%% assume the original signal is on [0,L].
%% assume the signal is on [0,1]. Frequencies are rescaled to \xi/L
xi = [(0:(n/2)) (((-n/2)+1):-1)]; xi = xi(:);
xhat = fft(x);
noctave = floor(log2(n)) - oct;
tfr = zeros(n, nvoice .* noctave);
kscale = 1;
tfrtic = zeros(1, nvoice .* noctave);
for jj = 1 : nvoice .* noctave
tfrtic(jj) = scale .* (2^(jj/nvoice));
end
for jo = 1:noctave % # of scales
for jv = 1:nvoice
qscale = scale .* (2^(jv/nvoice));
omega = xi ./ qscale ;
if strcmp(opts.motherwavelet, 'morse');
error('To use Morse wavelet, you have to download it...') ;
elseif strcmp(opts.motherwavelet, 'Cinfc');
tmp0 = (omega-opts.CENTER)./opts.FWHM;
tmp1 = (tmp0).^2-1;
windowq = exp( 1./tmp1 );
windowq( find( omega >= (opts.CENTER+opts.FWHM) ) ) = 0;
windowq( find( omega <= (opts.CENTER-opts.FWHM) ) ) = 0;
elseif strcmp(opts.motherwavelet, 'morlet')
windowq = 4*sqrt(pi)*exp(-4*(omega-0.69*pi).^2)-4.89098d-4*4*sqrt(pi)*exp(-4*omega.^2);
elseif strcmp(opts.motherwavelet, 'gaussian');
psihat = @(f) exp( -log(2)*( 2*(f-opts.CENTER)./opts.FWHM ).^2 );
windowq = psihat(omega);
elseif strcmp(opts.motherwavelet, 'meyer'); %% Meyer
windowq = zeros(size(omega));
int1 = find((omega>=5./8*0.69*pi)&(omega<0.69*pi));
int2 = find((omega>=0.69*pi)&(omega<7./4*0.69*pi));
windowq(int1) = sin(pi/2*meyeraux((omega(int1)-5./8*0.69*pi)/(3./8*0.69*pi)));
windowq(int2) = cos(pi/2*meyeraux((omega(int2)-0.69*pi)/(3./4*0.69*pi)));
elseif strcmp(opts.motherwavelet, 'BL3'); %% B-L 3
phihat = (2*pi)^(-0.5)*(sin(omega/4)./(omega/4)).^4; phihat(1) = (2*pi)^(-0.5);
aux1 = 151./315 + 397./840*cos(omega/2) + 1./21*cos(omega) + 1./2520*cos(3*omega/2);
phisharphat = phihat.*(aux1.^(-0.5));
aux2 = 151./315 - 397./840*cos(omega/2) + 1./21*cos(omega) - 1./2520*cos(3*omega/2);
aux3 = 151./315 + 397./840*cos(omega) + 1./21*cos(2*omega) + 1./2520*cos(3*omega);
msharphat = sin(omega/4).^4.*(aux2.^(0.5)).*(aux3.^(-0.5));
windowq = phisharphat.*msharphat.*exp(i*omega/2).*(omega>=0);
end
windowq = windowq ./ sqrt(qscale);
what = windowq .* xhat;
w = ifft(what);
tfr(:,kscale) = transpose(w);
kscale = kscale+1;
end
scale = scale .* 2;
end
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%% calculate the constant for reconstruction
%% TODO: calculate Rpsi for other mother wavelets
xi = [0.05:1/10000:10];
if strcmp(opts.motherwavelet, 'gaussian'); %% Gaussian (not really wavelet)
psihat = @(f) exp( -log(2)*( 2*(f-opts.CENTER)./opts.FWHM ).^2 );
windowq = psihat(xi);
Rpsi = sum(windowq./xi)/10000;
elseif strcmp(opts.motherwavelet, 'morlet');
windowq = 4*sqrt(pi)*exp(-4*(xi-0.69*pi).^2)-4.89098d-4*4*sqrt(pi)*exp(-4*xi.^2);
Rpsi = sum(windowq./xi)/10000;
elseif strcmp(opts.motherwavelet, 'Cinfc');
tmp0 = (xi - opts.CENTER)./opts.FWHM;
tmp1 = (tmp0).^2-1;
windowq = exp( 1./tmp1 );
windowq( find( xi >= (opts.CENTER+opts.FWHM) ) ) = 0;
windowq( find( xi <= (opts.CENTER-opts.FWHM) ) ) = 0;
Rpsi = sum(windowq./xi)/10000;
else
fprintf('Normalization is not implemented for Other mother wavelets, like BL3 and Meyer\n') ;
end
tfr = tfr ./ Rpsi;
end
實(shí)驗(yàn)結(jié)果
第一個(gè)為小波結(jié)果，第二個(gè)為同步壓縮小波變換結(jié)果。同步壓縮小波變換將小波系數(shù)擠壓到瞬時(shí)頻率處，所以它具有更加聚集的能量。小波變換的瑞利熵結(jié)果為13.1316，同步壓縮小波變換的瑞利熵結(jié)果為8.8856。
文章來源：http://www.zghlxwxcb.cn/news/detail-568602.html
[1] 何周杰. 同步變換理論、方法及其在工程信號(hào)分析中的應(yīng)用[D].上海交通大學(xué), 2020.
[2] ] I. Daubechies, J. Lu, H.-T. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2) (2011)243–261.文章來源地址http://www.zghlxwxcb.cn/news/detail-568602.html
到了這里，關(guān)于基于小波變換的同步壓縮變換原理和Matlab代碼的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
