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

%

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

This function is called by:

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

Generated on Tues 20-MAy-2008 10:40:08 by