• RMVANNUAL matlab remove annual cycle of a time series


    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

  • 相关阅读:
    Go并发编程实战 第2版 PDF (中文版带书签)
    DirectShow 应用开发过程
    Filter 原理
    DirectShow 常用函数总结
    COM 编程基础
    DirectShow 简介
    C++ 静态库与动态库以及在 Windows上 的创建、使用
    DirectShow 学习方法
    Qt 编译配置相关总结
    环境变量对于 VS 有什么用?
  • 原文地址:https://www.cnblogs.com/gisalameda/p/8288308.html
Copyright © 2020-2023  润新知