https://woodshole.er.usgs.gov/operations/sea-mat/sturges-html/rmvannual.html
% matlab file 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.
%
% ****************************
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