 clear;
 clc;
 % Prepare the independent variable (x-axis)of the synthesized data.
 % w_x is the corresponding spectral domain of x-space.
 % The synthesized data is a Gaussian profile modulating on a carrier with
 % a frequency of 20 and the phase is distorted by a cubic dispersion.
 x = linspace(0,1.0,256);
 % w_x=(0:255)'*1/1.0; 
 f_x=exp(-((x-0.5)/0.3).^2).*cos(2*pi*20*x).*exp(20*pi*sqrt(-1)*(x-0.5).^3);
 % 
 % phase2 contains the original phase profile 
 phase2=[];
 phase2=20*pi*(x-0.5).^3;
 phase2(:)=phase2(:)-phase2(128);
 figure(1)
 wfab=[];
 mwf=[];
 a=[1:1:60];
 % The modulated and phase distorted signal is analyzed with continuous
 % wavelet transform using Morlet wavelet cmor1-1 with a scaling spaning
 % from 1 to 40
 wfab=cwt(f_x,a,'cmor1-1','plot'); colormap jet;
 axis tight;
 mwf=abs(wfab);
 pwf=[];
 phase1=[];
 amplitude=[];
 pwf=angle(wfab);
 % The phase profile can be directly read from the phase pwf of the wavelet
 % tranform at the maximum of the modulus mwf along the scaling factor. 
 for m=1:256  % for each spectral components {m}
        maxwtm=mwf(1,m);
        k=1;
        for n=1:40
            if mwf(n,m)>maxwtm
                maxwtm=mwf(n,m);
                k=n;
            end
        end
        phase1(m)=pwf(k,m);
        amplitude(m)=mwf(k,m);
        ml(m)=k ;
 end
 phase1(:)=unwrap(phase1(:));
 phase1(:)=phase1(:)-phase1(128);
 yc10=[];
 xc10=[];
 % We remove the linear phase shift near the peak of the signal
 for jj=1:11
    yc10(jj)=phase1(round(128)+jj-6);
    xc10(jj)=x(round(128)+jj-6);
 end
 pc10=polyfit(xc10,yc10,1);
 phase1=phase1-polyval(pc10,x);
 phase1(:)=phase1(:)-phase1(128);
 % The red line is the modulus of the original signal
 % The blue line is the envelop of the signal retrieved from the Wavelet
 % transform.
 figure(2)
 subplot(211)
 axis([0,1,-1,1])
 plot(x,abs(f_x),'-r');
 hold on;
 scale=max(abs(f_x))/max(amplitude);
 % The phase profile is read out directly from the ridge of the wavelet
 % transform is shown by the blue symbols; while the original signal phase 
 % is shown by the red-colored curve 
 plot(x,scale.*amplitude,'-b'); 
 title('Magnitude of cwt along the modulus ridge');
 subplot(212)
 plot(x,phase1,'-bo'); 
 hold on;
 plot(x,phase2,'-r','LineWidth',2);
 title('Phase Profile along the modulus ridge');