% FFT example
clear variables
close all

set(0,'DefaultAxesFontSize',18)
set(0,'DefaultTextFontSize',18)
set(0,'DefaultLineLinewidth',1)


        title1= 'ENSO with  \Deltat=0.1[mon],Tmax=1000years'
        load ENSO_output.mat
        dt=delta_t/12;
        Fs=1/dt;
        Time=time;
Y=T;

    NTime=length(Time);
    
    clf
    set(gcf,'Position',[1700 1000 [600 1000]]);
    
    subplot(3,1,1)
    plot(Time,Y)
    xlabel 'Time [years]'
    title 'time series'
    ylabel('temperature anomaly [^\circC]')
    ylim([min(Y) max(Y)])
    text(2,1.4*max(Y),[title1])
    
    pause
    
    Yfft=fft(Y);
    ff=1/dt/NTime*(0:NTime-1);
    
    subplot(3,1,2)
    plot(ff,real(Yfft))
    hold on
    plot(ff,imag(Yfft))
    title('fft(Y)')
    xlabel('Frequency (#/year)')
    legend('Real','Imag')
    
    fpY=1/dt/NTime*(0:NTime/2);
    pY = (1/(Fs*NTime)) * abs(Yfft(1:NTime/2+1));
    pY(2:end-1) = 2*pY(2:end-1);
    
    pause 
    subplot(3,1,3)
    plot(fpY,pY(1:NTime/2+1))
    title('Power spectral density')
    xlabel('Frequency (#/year)')
    ylabel('deg^\circC/\Deltaf')
    
     pause
         yyaxis right
 semilogy(fpY,pY(1:NTime/2+1))

 xlim([0 2])
 

      Nsegments=5;
      for ii=1:Nsegments
          NTime1=NTime/Nsegments;
          Y1=Y((ii-1)*NTime1+1:(ii)*NTime1);
    %      Y1=Y1-mean(Y1);
    Yfft=fft(Y1);
    ff=1/dt/NTime1*(0:NTime1-1);
    
    pY = (1/(Fs*NTime1)) * abs(Yfft(1:NTime1/2+1));
    pY(2:end-1) = 2*pY(2:end-1);
    pY_all(ii,:)=pY;
      end
          fpY=1/dt/NTime1*(0:NTime1/2);

      pY=mean(pY_all,1);
    
    subplot(3,1,3)
    hold on
    semilogy(fpY,pY(1:NTime1/2+1),'k','Linewidth',2)
    
    
    pause
    
    figure
    
    subplot(3,1,1)
       findpeaks(Y,Time,'MinPeakProminence', 3);

    [pks,locs] =findpeaks(Y,'MinPeakProminence', 3);

    enso_int=nearest(diff(locs)/Fs);
    
    subplot(3,1,2)
    histogram(enso_int)
    
    subplot(3,1,3)
    
        time_int=-1:dt:1;
        Y_event=zeros(1,2*Fs+1);
    for ii=3:length(locs)
        itime_int=locs(ii)-1*Fs:locs(ii)+1*Fs;
        Y_event=Y_event+Y(itime_int);
        plot(time_int,Y(itime_int))
        hold on
    end
    Y_event=Y_event/(length(locs)-1);
 %   plot(time_int,Y_event,'k','LineWidth',5)
    title Composites
   
    
    
