Deformation estimations on a synthetic signal
Copyright (C) 2017 Adrien MEYNARD
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.
Author: Adrien MEYNARD Email: adrien.meynard@univ-amu.fr Created: 2017-12-19
Contents
Load signal
clear all; close all; clc; warning off; addpath('cwt'); addpath('deform_estimation'); load('signals/toy_sig_paper'); T = length(y);
Joint estimation
Dt = 200; % temporal subsampling for the deformation estimation dgamma0 = ones(1,T); % gamma'(t) initialization a0 = ones(1,T); % a(t) initialization wav_typ = 0; % wavelet type (cf. cwt.m) wav_paramWP = 25; % corresponding parameter for warping estimation wav_param = 500; % corresponding parameter for spectrum and AM estimations NbScales = 106; scalesAM = 2.^(linspace(-1,3.5,NbScales)); subrate = 7; % subsampling step for the scales to ensure the covariance invertibility scalesWP = scalesAM(1:subrate:end); stopWP = 2e-2; % minimal gap between two steps in the gradient itWP = 6; % number of gradient iterations r = 1e-3; % regularization parameter Nf = 2500; % number of frequencies for spectrum estimation NbScalesS = 212; scalesS = 2.^(linspace(-1,3.5,NbScalesS)); % for spectrum estimation Nit = 10; % maximal number of iterations in the joint estimation stop_crit = 1e-3; % relative update threshold paramWAV = {wav_typ,wav_param,wav_paramWP}; paramAM = {'AM',scalesAM,r}; % AM (model without noise) paramWP = {scalesWP,itWP,stopWP}; paramS = {scalesS,Nf}; tic; [aML, dgammaML, Sxest, evol_crit] = estim_altern(y,Dt,dgamma0,a0,paramWAV,paramWP,paramAM,paramS,stop_crit,Nit); toc;
Iteration 1 Relative update WP: Inf % Relative update AM: 71.71 % Iteration 2 Relative update WP: 34.37 % Relative update AM: 33.32 % Iteration 3 Relative update WP: 35.52 % Relative update AM: 105.53 % Iteration 4 Relative update WP: 71.20 % Relative update AM: 72.27 % Iteration 5 Relative update WP: 1.75 % Relative update AM: 11.91 % Iteration 6 Relative update WP: 0.11 % Relative update AM: 0.28 % Iteration 7 Relative update WP: 0.01 % Relative update AM: 0.07 % Elapsed time is 69.789505 seconds.
Deformations estimation
addpath('analysis') Wy = cwt(y,scalesAM,0,wav_param); thetaWP_B = baselinewarpest(Wy,scalesAM); thetaAM_B = baselineAMest(Wy); thetaWPV = log2(dgamma(1:Dt:end)).'; thetaAMV = (a(1:Dt:end).^2).'; bWP = crlbWP(thetaWPV,thetaAMV,Sx+eps,scalesWP,wav_typ,wav_paramWP); WPinf = thetaWPV - 3*sqrt(bWP); WPsup = thetaWPV + 3*sqrt(bWP); bAM = crlbAM(thetaAMV,scalesAM); AMinf = thetaAMV - 3*sqrt(bAM); AMsup = thetaAMV + 3*sqrt(bAM); t = linspace(0,(T-1)/Fs,T); th = t(1:Dt:end); figure; plot(t,thetaAM_B,'g',t,aML.^2,'r',t,a.^2,'b-.',th,AMinf,'b--',th,AMsup,'b--','linewidth',2); grid on; title('Amplitude modulation estimation'); legend({'Baseline estimator','Algorithm estimator','$\theta_{n,1}$','$\theta_{n,1} \pm 3\sqrt{CRLB(\theta_{n,1})}$'},'Interpreter','latex'); axis([0 t(end) 0 3]); xlabel('Time (s)'); ylabel(''); figure; plot(t,thetaWP_B,'g',t,log2(dgammaML),'r',t,log2(dgamma),'b-.',th,WPinf,'b--',th,WPsup,'b--','linewidth',2); grid on; axis tight; title('Warping estimation'); legend({'Baseline estimator','Algorithm estimator','$\theta_{n,2}$','$\theta_{n,2}\pm 3\sqrt{CRLB(\theta_{n,2})}$'},'interpreter','latex'); xlabel('Time (s)'); ylabel('');
![](estim_WPAM_toysig_01.png)
![](estim_WPAM_toysig_02.png)
Spectrum estimation
freq = linspace(0,Fs,2*Nf - 1); figure; plot(freq,Sx,'-.b',freq,Sxest,'r','linewidth',2); xlabel('Frequency (Hz)');ylabel('Spectrum'); grid on; axis tight; legend('True spectrum','Estimated spectrum'); V=axis; axis([300 1600 0 V(4)]);
![](estim_WPAM_toysig_03.png)
Stopping criterion
Nit = size(evol_crit,1); figure; semilogy(evol_crit,'linewidth',2); hold on; semilogy(1:Nit,0.001*ones(Nit,1),'k--','linewidth',2); % threshold text(1.5,1.5e-3,'\fontsize{30}T = 0.1 %'); legend('Time warping','Amplitude modulation','Spectrum'); axis tight; grid on; xlabel('Iteration'); ylabel('Relative update');
![](estim_WPAM_toysig_04.png)
MSE
errAM = (aML(:).^2 - a.^2).^2; errAM = errAM(Dt:end); % First sample of thetaAM is rejected because of boundaries effect mseAM = mean(errAM.^2); mseAM0 = mean((thetaAM_B(:) - a.^2).^2); mseWP = mean((log2(dgammaML(:)) - log2(dgamma)).^2); mseWP0 = mean((thetaWP_B(:) - log2(dgamma)).^2); fprintf('\nMean square errors: \n\n AM estimation:\n Baseline method: %.2e \n Algorithm estimation: %.2e \n\n Time Warping estimation:\n Baseline estimation: %.2e \n Algorithmestimation: %.2e\n',[mseAM0,mseAM,mseWP0,mseWP])
Mean square errors: AM estimation: Baseline method: 2.01e-01 Algorithm estimation: 7.01e-02 Time Warping estimation: Baseline estimation: 2.32e-02 Algorithmestimation: 4.91e-04