% %   This file is part of the TFPM toolbox v0.5 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Playing round with real data... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ORIGINAL DATAload('/users/mjachan/tex/prj/TF/yearly')Sun= yearly(61, 1:2);Sun= [Sun; yearly(1:60, 1: 2)];Sun= [Sun; yearly(1:60, 3: 4)];Sun= [Sun; yearly(1:60, 5: 6)];Sun= [Sun; yearly(1:60, 7: 8)];Sun= [Sun; yearly(1:60, 9:10)];Sun= [Sun; yearly(61, 3:4)];Sun= [Sun; yearly(61, 5:6)];Sun= [Sun; yearly(61, 7:8)];N0= size(Sun, 1);sun0= Sun(:, 2);% TIME SERIES PLOTfigure(1);clf;hold onplot(Sun(:, 1), sun0, 'Color', 'k', 'Linewidth', 2, 'Linestyle', '-')axis([1700, 2002, 0, 200])%ticksoff;% MISSING DATAmissing= sun0==0;Sun(missing)% VERY SHORT TIME AVERAGINGNwin= 1;window= hanning(2*Nwin+1);window= window/sum(window);sun0short= conv(window, sun0);sun0short= sun0short(Nwin+1:Nwin+N0);% RESTORE MISSING DATAsun1= sun0+sun0short.*missing;% MEDIUM-SHORT TIME AVERAGINGNwin= 15;window= hanning(2*Nwin+1);window= window/sum(window);sun1short= conv(window, sun1);sun1short= sun1short(Nwin+1:Nwin+N0);plot(Sun(:, 1), sun1short, 'Color', 'k', 'Linewidth', 1, 'Linestyle', '--')set(gca, 'YTick', [40 80 120 160 200])boxprint -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_orig% HISTOGRAM PLOTfigure(2);clf;hold onHsun1= hist(sun1, 10:20:190);barh(10:20:190, Hsun1, 1)axis([0, 125, 0, 200])colormap(gray/2+.5)% FIT AN EXPONENTIAL DENSITYMAXsun= 200%!!!lambda= 1/mean(sun1);xi= 1:MAXsun;fxi= N0*MAXsun/10*lambda*exp(-lambda.*xi);plot(fxi, xi, 'Color', 'k', 'Linewidth', 3, 'Linestyle', '-')set(gca, 'XTick', [0 20 40 60 80 100 120])set(gca, 'YTick', [10 30 50 70 90 110 130 150 170 190])box%ticksoff;print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_orig_hist% SUBTRACT SHORT-TIME MEANsun2= sun1-sun1short;% TIME SERIES PLOTfigure(3);clf;hold onplot(sun2, 'Color', 'k', 'Linewidth', 2, 'Linestyle', '-')axis([0, N0-1, -120, 120])set(gca, 'YTick', [-120 -80 -40  0 40 80 120])%ticksoff;boxprint -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_gaus% HISTOGRAM PLOTfigure(4);clf;hold onHsun2= hist(sun2, -108:24:108);barh(-108:24:108, Hsun2, 1)axis([0, 90, -120, 120])colormap(gray/2+.5)% FIT A GAUSSIAN DENSITYsi= std(sun2);MAXsun= 120;%!!!!xi= -MAXsun:1:MAXsun;fxi= 2*MAXsun*N0/10*exp(-xi.^2/2/si^2)/sqrt(2*pi)/si;plot(fxi, xi, 'Color', 'k', 'Linewidth', 3, 'Linestyle', '-')set(gca, 'XTick', [0 20 40 60 80])set(gca, 'YTick', [-108   -84   -60   -36   -12    12    36    60    84   108])box%ticksoff;print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_gaus_hist% CORRELATIONSx= sun2;qB= xcorr(x, x, 'biased');qU= xcorr(x, x, 'unbiased');qb= qB(N0:end);qu= qU(N0:end);qc= real(ifft(abs(fft(x)).^2))/N0;figure(5);clf;hold onplot(qb, 'Color', 'k', 'Linewidth', 1, 'Linestyle', '-')plot(qu, 'Color', 'k', 'Linewidth', 2, 'Linestyle', ':')plot(qc, 'Color', 'k', 'Linewidth', 1, 'Linestyle', '--')axis([0 N0-1 -1500 1500])box%ticksoff;print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_corr% PERIODOGRAMPER= real(fft([qB(N0:end); qB(1:N0-1)]));COR= real(fft([qU(N0:end); qU(1:N0-1)]));figure(6);clf;hold onplot(log10(PER(1:N0)), 'Color', 'k', 'Linewidth', 2, 'Linestyle', '-')%plot(COR(1:N0), 'Color', 'g', 'Linewidth', 1, 'Linestyle', '-')% WINDOWED PERIODOGRAMNwin= 99;% ODD!!!!!!!!window= [zeros(N0-(Nwin+1)/2, 1); hamming(Nwin); zeros(N0-(Nwin+1)/2, 1)];rB= qB.*window;rU= qU.*window;PERw= real(fft([rB(N0:end); rB(1:N0-1)]));CORw= real(fft([rU(N0:end); rU(1:N0-1)]));plot(log10(PERw(1:N0)), 'Color', 'r', 'Linewidth', 2, 'Linestyle', '-')Nwin= 49;% ODD!!!!!!!!window= [zeros(N0-(Nwin+1)/2, 1); hamming(Nwin); zeros(N0-(Nwin+1)/2, 1)];rB= qB.*window;rU= qU.*window;PERw= real(fft([rB(N0:end); rB(1:N0-1)]));CORw= real(fft([rU(N0:end); rU(1:N0-1)]));plot(log10(PERw(1:N0)), 'Color', 'g', 'Linewidth', 2, 'Linestyle', '-')axis([0 N0-1 -1 5])box%ticksoff;print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_peri%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBSAMPLING BY 2sun= sun3(1:2:end);N0= length(sun);% FIT INTO 256 FRAME WITH NOISEN= 256;SNR= 10;%dBx= [zeros((N-N0)/2, 1); sun; zeros((N-N0+1)/2, 1)]+0.1*randn(N, 1);plot(x)alpha= 1/2;%%%% NONPARAMETRIC ANALYSIS% PERIODOGRAMfigure(7);clf;hold onperiodogram(sun, [], 512, 'twosided')pyulear(sun, 10, [], 'twosided')% SPECTROGRAMNwin= 32;figure(8);clfspecgram(sun, Nwin, [], [], Nwin-1)%%%% PARAMETRIC ANALYSISMARMAX= 8;LARMAX= 3;MAR= 5;LAR= 1;SIG= zeros(MARMAX, LARMAX+1);for MAR= 1:MARMAX   for LAR= 0:LARMAX      Psi= conj(tf_multiwin(N, MAR, LAR, 6, 4, 1));      psi= Psi(N/2-3*LAR+1:N/2+3*LAR+1, N/2-MAR+1:N/2+MAR+1);      Qxx= corr_est(x, x, MAR, alpha);      Axx= fft(Qxx);      Axx= [Axx(N/2+1:N, :); Axx(1:N/2, :)];      axx= Axx(N/2+1-3*LAR:N/2+1+3*LAR, :);      [AMLw, B0Lw]= tfar_est_tfywu(axx.*psi, N);      Amlw= AMLw(:, :, end);      B0lw= B0Lw(:, :, end);      Pxxw= tfarma_wvsp(Amlw, B0lw, N, alpha);      e= real(tfarma_inv(x, Amlw, B0lw));      SIG(MAR, LAR+1)= var(e);      figure(1);clf;plot(e, 'Color', 'k', 'Linewidth', 1, 'Linestyle', '-')      if(LAR<2)	 axis([0 255 -6 6])         text(150, 5, sprintf('(%d, %d); %0.5g', MAR, LAR, var(e)))      else	 axis([0 255 -50 50])         text(150, 40, sprintf('(%d, %d); %0.5g', MAR, LAR, var(e)))      end      resname= sprintf('sunInn%02d%02d', MAR, LAR);      drawnow      print('-depsc2', strcat('/users/mjachan/tex/prj/TF/figures/', resname));            nk_plot(Pxxw, 2);      colormap(flipud(gray/2+.5))      resname= sprintf('sunWxx%02d%02d', MAR, LAR);      text(150, 30, sprintf('(%d, %d)', MAR, LAR))      drawnow      print('-depsc2', strcat('/users/mjachan/tex/prj/TF/figures/', resname));%      [AMLu, B0Lu]= tfar_est_tfywu(axx, N);%      Amlu= AMLu(:, :, end);%      B0lu= B0Lu(:, :, end);%      Pxxu= tfarma_wvsp(Amlu, B0lu, N, alpha);%      nk_plot(Pxxu, 3);%      resname= sprintf('sunPxx%02d%02d', MAR, LAR);%      title(resname)%      drawnow%      print('-depsc2', strcat('/users/mjachan/tex/prj/TF/figures/', resname));   endend[MAR, LAR]=find(SIG==min(min(SIG)));LAR= LAR-1;Psi= conj(tf_multiwin(N, MAR, LAR, 6, 4, 1));psi= Psi(N/2-3*LAR+1:N/2+3*LAR+1, N/2-MAR+1:N/2+MAR+1);Qxx= corr_est(x, x, MAR, alpha);Axx= fft(Qxx);Axx= [Axx(N/2+1:N, :); Axx(1:N/2, :)];axx= Axx(N/2+1-3*LAR:N/2+1+3*LAR, :);[AMLw, B0Lw]= tfar_est_tfywu(axx.*psi, N);Amlw= AMLw(:, :, end);B0lw= B0Lw(:, :, end);P= tfarma_wvsp(Amlw, B0lw, N, alpha);H= real(tfarma_impr(Amlw, B0lw, N, alpha));S= tfarma_sprf(Amlw, B0lw, N, alpha);R= real(tfarma_corr(Amlw, B0lw, N, alpha));clear;tfpm;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%