【数字信号处理】小波变换的时间重分配多重同步压缩变换【含Matlab源码 2460期】.zip

function [Ws,Wx,tau2,as,fs] = WMTSST_Z(s,gamma,mywav,totalscal,num) % function Wsst : computes the time-reassigned multisynchrosqueezing transform of signal s. % Implements algorithm of Brevdo, Daubechies, Wu et al. % % Inputs % s : input signal % gamma : threshold % mywav : string coding the mother wavelet % totalscal : number of discretized scales % % Outputs : % Ws : the time-reassigned multisynchrosqueezing transform % Wx : the wavelet transform % tau : group delay estimator % as : scales vector % fs : frequency vector if nargin<5 gamma = 0.01; mywav = 'cmor2-1'; totalscal=256; num=10; end % Computing scales (0-fs/2) wcf=centfrq(mywav); cparam=2*wcf*totalscal; a=totalscal:-1:1; as=cparam./a; [q,na]=size(a); x = s(:); n = length(s); %% Wavelet transform and reassignment operators Wx = zeros(na, n); % CWT tau = zeros(na, n); % Horizontal extimator x = x(:).'; xh = fft(x); xi = (0:n-1)/n; % Filter definition if strncmp(mywav,'gmor',4) [v1 v2] = regexp(mywav,'[0-9.]*-[0-9.]'); beta = str2num(mywav(v1:v2-2)); gam = str2num(mywav(v2:end)); filt = @(a) gmor(beta,gam,a*xi); func = @(x) 2*(exp(1)*gam/beta)^(beta/gam) * x.^beta .* exp(-x.^gam); elseif strncmp(mywav,'cmor',4) [v1 v2] = regexp(mywav,'[0-9.]*-[0-9.]'); Fb = str2num(mywav(v1:v2-2)); Fc= str2num(mywav(v2:end)); filt = @(a) cmor(Fb,Fc,a*xi); func = @(x) sqrt(Fb) * exp(-Fb^2*pi*(x-Fc).^2); elseif strncmp(mywav,'bump',4) [v1 v2] = regexp(mywav,'[0-9.]*-[0-9.]'); mu = str2num(mywav(v2:end)); sigma = str2num(mywav(v1:v2-2)); filt = @(a) bump(mu,sigma,a*xi); func = @(x) exp(1-1./(1-((x-mu)./sigma).^2)); end % for each octave for ai = 1:length(as) a = as(ai); [psih dfilt] = filt(a); wxtmp = ifft(conj(psih) .* xh); %WT Wx(ai,:) = wxtmp; tau(ai, :) = a*real(ifft(conj(dfilt).*xh) ./ wxtmp/2/1i/pi); end Wx = Wx(:, 1:n); tau = tau(:, 1:n); [M,H]=size(tau); tau2= zeros(M,H); for j=1:M tau(j,:)=(1:n)+tau(j,:); end [neta,nb]=size(Wx); if num>1 for kk=1:num-1 for b=1:nb for eta=1:neta k2 = round(tau(eta,b)); if k2>=1 && k2<=nb tau2(eta,b)=tau(eta,k2); end end end tau=tau2; end else tau2=tau; end tau2=round(round(tau2*2)/2); tau2(abs(Wx)<gamma) = NaN; fs = Fc./as; Ws = zeros(size(Wx)); Wx2 = zeros(size(Wx)); [mm,nn]=size(Wx); for eta1=1:mm%frequency Wx2(eta1,:)=Wx(eta1,:).*exp(-1j* 2 * pi*fs(eta1)*(1:nn)); end for b=1:n for ai=1:length(as) k=round(tau2(ai,b)); if k>0 && k<=n Ws(ai,k) = Ws(ai, k) + Wx2(ai, b); end end end end