RMVANNUAL

% 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.

%

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

% 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.

%

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

0001 %  matlab file  RMVANNUAL
0002 %
0003 %  Purpose is to remove the annual cycle of a time series -- carefully.
0004 %      The method is to compute the  FFT,  replace the power at the annual
0005 %      frequency by averaging the values on each side, rather than removing
0006 %      it entirely.  Then the fft is inverted to return to the time domain.
0007 %      The same thing is done in phase
0008 %
0009 %      The program must be run again for each harmonic; it seems
0010 %      appropriate to avoid mucking up the time series if there is little power
0011 %      in the vicinity of a specific harmonic term, altho' it is a bit of a
0012 %      bother (but trivial) to have to run it several times.  
0013 %
0014 %      The program works equally well for any specific Fourier term of
0015 %      course
0016 %
0017 %   5 May 2008   WS
0018 %      
0019 %      version 1.0  beta
0020 %
0021 %    **************************
0022 %     PROGRAM REQUIRES THAT THE INPUT FILE BE NAMED  " X "
0023 %     ---  in lower case ---
0024 %     and calls an external routine  DOPGRAM 
0025 %     N.B.  the program does not check to be sure the length of the data
0026 %     vector is an exact multiple of the annual frequency; user's job
0027 %    
0028 %     The output vector is  "  X2 " ; input x is detrended.
0029 %
0030 %     ****************************
0031 x = detrend(x);
0032 %  clear the graphics window just to be safe
0033 clf
0034 %
0035 %   First we compute the Fourier transform and look to see if there is a
0036 %   large amount of power at the frequencies we are concerned about
0037 %    use a routine called dopgram; it only shows power in the raw
0038 %    periodogram
0039 %
0040 xin = x;
0041 dopgram
0042 loglog(rawpgram)
0043 xlabel (' raw fft --hit any key to continue')
0044 pause
0045 %
0046 %  compute the Fourier transform and plot power and phase
0047 pin = fft(x);
0048 pinsave=pin;
0049 subplot(211),plot(abs(pin)), title('  Power')
0050 xlabel(' sine terms come first, then cosine terms ')
0051 subplot(212),plot(angle(pin)/pi),title ('  Phase, hit any key to continue'), pause
0052 %
0053 l = length(pin)     % just to show the length of the signal
0054 nout = input ('   Which frequency component do we want suppressed?     ')
0055 % smooth over the adjacent bands
0056 nn=nout+1 ; %  the first term is the mean value 
0057 pin(nn) = ( pin(nn-1) + pin(nn+1))/2;   % that gets the "a" term - sine
0058 n2=l-nout+1 ;
0059 pin(n2) = ( pin(n2-1) + pin(n2+1))/2;   %   and the "b" term - cosine
0060 % now look at it to be safe
0061 %  but only plot a limited part of Freq space for high resolution
0062 clf
0063 plot(abs(pin(1:3*nn)),'-x'), title ( '    Power after the freq point has been clipped  ')
0064 hold on, plot(abs(pinsave(1:3*nn)),'r')
0065 xlabel(' original is in red, hit any key to continue')
0066 pause
0067clear x2, hold off
0068 x2 = real(ifft(pin));
0069 %
0070 %     Now that we've done it, let's look at it.
0071 %     first, look at the fft
0072 pgsave=rawpgram;
0073 xin=x2;
0074 dopgram
0075 loglog(pgsave),hold on, loglog(rawpgram,'r')
0076 xlabel( ' the clipped version is in red, hit any key to continue')
0077 pause
0078 clf
0079 % now look at the time series itself
0080 hold off, plot(x), hold, plot(x2, 'g')
0081 title ( 'Signal before & after (G) the freq point was  clipped  ')
0082 xlabel (' hit any key to continue' ),pause
0083 hold off, clf