tfpm_sunspot.m

用于模拟时变非平稳的ARMA过程

% 
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] 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 DATA
load('/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 PLOT
figure(1);clf;hold on
plot(Sun(:, 1), sun0, 'Color', 'k', 'Linewidth', 2, 'Linestyle', '-')
axis([1700, 2002, 0, 200])
%ticksoff;

% MISSING DATA
missing= sun0==0;
Sun(missing)

% VERY SHORT TIME AVERAGING
Nwin= 1;
window= hanning(2*Nwin+1);
window= window/sum(window);
sun0short= conv(window, sun0);
sun0short= sun0short(Nwin+1:Nwin+N0);

% RESTORE MISSING DATA
sun1= sun0+sun0short.*missing;

% MEDIUM-SHORT TIME AVERAGING
Nwin= 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])
box
print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_orig

% HISTOGRAM PLOT
figure(2);clf;hold on
Hsun1= hist(sun1, 10:20:190);
barh(10:20:190, Hsun1, 1)
axis([0, 125, 0, 200])
colormap(gray/2+.5)

% FIT AN EXPONENTIAL DENSITY
MAXsun= 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 MEAN
sun2= sun1-sun1short;

% TIME SERIES PLOT
figure(3);clf;hold on
plot(sun2, 'Color', 'k', 'Linewidth', 2, 'Linestyle', '-')
axis([0, N0-1, -120, 120])
set(gca, 'YTick', [-120 -80 -40  0 40 80 120])
%ticksoff;
box
print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/sun_gaus

% HISTOGRAM PLOT
figure(4);clf;hold on
Hsun2= hist(sun2, -108:24:108);
barh(-108:24:108, Hsun2, 1)
axis([0, 90, -120, 120])
colormap(gray/2+.5)

% FIT A GAUSSIAN DENSITY
si= 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

% CORRELATIONS
x= 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 on
plot(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

% PERIODOGRAM
PER= real(fft([qB(N0:end); qB(1:N0-1)]));
COR= real(fft([qU(N0:end); qU(1:N0-1)]));

figure(6);clf;hold on
plot(log10(PER(1:N0)), 'Color', 'k', 'Linewidth', 2, 'Linestyle', '-')
%plot(COR(1:N0), 'Color', 'g', 'Linewidth', 1, 'Linestyle', '-')

% WINDOWED PERIODOGRAM
Nwin= 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 2
sun= sun3(1:2:end);
N0= length(sun);

% FIT INTO 256 FRAME WITH NOISE
N= 256;
SNR= 10;%dB
x= [zeros((N-N0)/2, 1); sun; zeros((N-N0+1)/2, 1)]+0.1*randn(N, 1);
plot(x)
alpha= 1/2;

%%%% NONPARAMETRIC ANALYSIS

% PERIODOGRAM
figure(7);clf;hold on
periodogram(sun, [], 512, 'twosided')

pyulear(sun, 10, [], 'twosided')


% SPECTROGRAM
Nwin= 32;
figure(8);clf
specgram(sun, Nwin, [], [], Nwin-1)


%%%% PARAMETRIC ANALYSIS
MARMAX= 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));
   end
end
[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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%