RMVANNUAL

 %  matlab file  RMVANNUAL

 %http://woodshole.er.usgs.gov/operations/sea-mat/sturges-html/index.html

 %  Purpose is to remove the annual cycle of a time series -- carefully.

 %      The method is to compute the  FFT,  replace the power at the annual

 %      frequency by averaging the values on each side, rather than removing

 %      it entirely.  Then the fft is inverted to return to the time domain.

 %      The same thing is done in phase

 %

 %      The program must be run again for each harmonic; it seems

 %      appropriate to avoid mucking up the time series if there is little power

 %      in the vicinity of a specific harmonic term, altho' it is a bit of a

 %      bother (but trivial) to have to run it several times.  

 %

 %      The program works equally well for any specific Fourier term of

 %      course

 %

 %   5 May 2008   WS

 %      

 %      version 1.0  beta

 %

 %    **************************

 %     PROGRAM REQUIRES THAT THE INPUT FILE BE NAMED  " X "

 %     ---  in lower case ---

 %     and calls an external routine  DOPGRAM 

 %     N.B.  the program does not check to be sure the length of the data

 %     vector is an exact multiple of the annual frequency; user's job

 %    

 %     The output vector is  "  X2 " ; input x is detrended.

 %

 %     ****************************

 x=sw;

 x = detrend(x);

 %  clear the graphics window just to be safe

 clf

 %

 %   First we compute the Fourier transform and look to see if there is a

 %   large amount of power at the frequencies we are concerned about

 %    use a routine called dopgram; it only shows power in the raw

 %    periodogram

 %

 xin = x;

 dopgram

 loglog(rawpgram)

 xlabel (' raw fft --hit any key to continue')

 pause

 %

 %  compute the Fourier transform and plot power and phase

 pin = fft(x);

 pinsave=pin;

 subplot(211),plot(abs(pin)), title('  Power')

 xlabel(' sine terms come first, then cosine terms ')

 subplot(212),plot(angle(pin)/pi),title ('  Phase, hit any key to continue'), pause

 %

 l = length(pin)     % just to show the length of the signal

 nout = input ('   Which frequency component do we want suppressed?     ')

 % smooth over the adjacent bands

 nn=nout+1 ; %  the first term is the mean value 

 pin(nn) = ( pin(nn-1) + pin(nn+1))/2;   % that gets the "a" term - sine

 n2=l-nout+1 ;

 pin(n2) = ( pin(n2-1) + pin(n2+1))/2;   %   and the "b" term - cosine

 % now look at it to be safe

 %  but only plot a limited part of Freq space for high resolution

 clf

 plot(abs(pin(1:3*nn)),'-x'), title ( '    Power after the freq point has been clipped  ')

 hold on, plot(abs(pinsave(1:3*nn)),'r')

 xlabel(' original is in red, hit any key to continue')

 pause

clear x2, hold off

 x2 = real(ifft(pin));

 %

 %     Now that we've done it, let's look at it.

 %     first, look at the fft

 pgsave=rawpgram;

 xin=x2;

 dopgram

 loglog(pgsave),hold on, loglog(rawpgram,'r')

 xlabel( ' the clipped version is in red, hit any key to continue')

 pause

 clf

 % now look at the time series itself

 hold off, plot(x), hold, plot(x2, 'g')

 title ( 'Signal before & after (G) the freq point was  clipped  ')

 xlabel (' hit any key to continue' ),pause

 hold off, clf