ipf.m

这是一个时间序列的高通滤波matlab函数！

function ipf(x,y)
% ipf(x,y)
% Keyboard-operated Interactive Peak Fitter for data in arguments x,y.
% This version does not use sliders; it uses the keyboard commands only.
% See http://www.wam.umd.edu/~toh/spectrum/InteractivePeakFitter.htm
% T. C. O'Haver ([e-mail]).  Version 3.5: April 2008. Model is computed 
% at 100 points regardless of number of data points in version 3.2 and
% later). 
% Example: x=[0:.005:1];y=humps(x).^3;ipf(x,y)
% For large data sets, to view only a portion of the data over a specified
% x-axis range, you can type n1=val2ind(x,x1);n2=val2ind(x,x2);ipf(x(n1:n2),y(n1:n2))
% where x1 and x2 are the end limits of the x-axis values displayed.
%
% Keyboard Controls:
% Pan signal left and right:  Coarse pan: < and >   
%                             Fine pan: left and right cursor arrow keys
% Zoom in and out:            Coarse zoom: / and '   
%                             Fine zoom: up and down cursor arrow keys
% Select # of peaks:          Number keys 1-5
% Select peak shape:          g Gaussian
%                             l Lorentzian
%                             o Logistic
%                             p Pearson (use a,z keys to adjust shape)
%                             e exponentially-broadened Gaussian 
%                                 (use a,z keys to adjust broadening)
% Fit:                        f
% Toggle autozero off/on      t   (Newly added as of version 3.5)
% Baseline:                   b, then click left and right baseline
% Custom start position:      c, then click on each peak position
% Adjust 'extra' up or down:  a,z
% Print parameters & results: q
% Print fit results only:     r
% Print out data segment:     d  (Newly added as of version 3.3)

global X
global Y
global xo
global dx
global NumPeaks
global Shape
global MeanFitError
global PEAKHEIGHTS
global FitResults
global start
global extra
global delta
global AUTOZERO
close
format short g
format compact
warning off all

% Assign arguments to internal global variables
X=x;
Y=y;

% Adjust X and Y vector shape to 1 x n (rather than n x 1)
X=reshape(X,1,length(X));
Y=reshape(Y,1,length(Y));

% X=[1:length(Y)]; % Use this only to create an x vector if needed

% Remove excess offset from data
Y=Y-min(Y);

% Initial values of parameters
NumPeaks=1; % Initial Number of peaks in model
Shape=1; % Initial Shape of the peaks (1=Gaussian, 2=Lorentzian, etc)
extra=.5;
delta=0; % Initial Random change in initial start guesses for each re-fit
xo=length(Y)/2; % Initial Pan setting
dx=length(Y)/4; % Initial Zoom setting
start=length(Y)/2;
PEAKHEIGHTS=zeros(NumPeaks,1);
AUTOZERO=1; % Change to 0 to disable autozero operation

% Plot the signal and its fitted components and residuals
RedrawSignal(X,Y,xo,dx,NumPeaks);
h=figure(1);
h2=gca;
% Attaches KeyPress test function to the figure.
set(gcf,'KeyPressFcn',@ReadKey)
uicontrol('Style','text')
% end of outer function
% ----------------------------SUBFUNCTIONS--------------------------------
function ReadKey(obj,eventdata)
% Interprets key presses from the Figure window.
global X
global Y
global xx
global yy
global xo
global dx
global start
global FitResults
global NumPeaks
global Shape
global delta
global PEAKHEIGHTS
global AUTOZERO
global extra
global MeanFitError
% When a key is pressed, interprets key and calls corresponding function.
% Note: If you don't like my key assignments, you can change the numbers
% in the case statements here to re-assign that function to any other key.
key=get(gcf,'CurrentCharacter');
if ischar(key),
  switch double(key),
    case 28
        % Pans one point down when left arrow pressed.
          xo=xo-1;
          if xo<1,xo=1;,end
          [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
         % h2=gca;
    case 29
        % Pans one point up when right arrow pressed.
        ly=length(Y);
        xo=xo+1;
        if xo>ly,xo=ly;,end
        [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
    case 44
        % Pans 2% down when < key pressed.
        ly=length(Y);
        xo=xo-ly/50;
        if xo<1,xo=1;,end
        [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
    case 46
        % Pans 2% up when > key pressed.
        ly=length(Y);
        xo=xo+ly/50;
        if xo>ly,xo=ly;,end
        [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
    case 31
        % Zooms one point up when up arrow pressed.
        dx=dx+2;
        [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
    case 30
        % Zooms one point down when down arrow pressed.
        dx=dx-2;
        [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
    case 47
        % Zooms 2% up when / pressed.
        ly=length(Y);
        dx=dx+ly/50;
        [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
    case 39
        % Zooms 2% down when ' pressed.
        ly=length(Y);
        dx=dx-ly/50;
        [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);      
    case 102
        % When 'f' key is pressed, computes and plots fit.
        startnow=start;
        delta=(max(xx)-min(xx))/100;
          for k=1:2:2*NumPeaks,
            startnow(k)=start(k)+randn*delta;
          end
        [FitResults,MeanFitError]=FitAndPlot(xx,yy,NumPeaks,Shape,delta,startnow,extra); 
    case 99
        % When 'c' key is pressed, user clicks graph to enter start positons, 
        % then fit is computed and graph re-drawn.q
        % Acquire first-guess peak positions from user mouse pointer
        subplot(2,1,1);xlabel('Click on the estimated positions of each proposed component peak.')
        [clickX,clickY] = GINPUT(NumPeaks);
        % Create a "start" vector using these peak positions, with peak widths equal 
        % to 1/10 of the zoom region.
        n=max(xx)-min(xx);
        width=n/(5*NumPeaks);
        start=[];
        for k=1:NumPeaks,
            start=[start clickX(k) width];
        end
        [FitResults,MeanFitError]=FitAndPlot(xx,yy,NumPeaks,Shape,delta,start,extra);    
    case 98
        % When 'b' key is pressed, user clicks graph before and after peaks
        % to enter background points, then fit is computed and graph re-drawn.
          subplot(2,1,1);xlabel('Click on the baseline to the LEFT the peak(s).')
          [X1,Y1] = GINPUT(1);
          subplot(2,1,1);xlabel('Now click on the baseline to the RIGHT the peak(s).')
          [X2,Y2] = GINPUT(1);
          n=length(xx);
          %  Create "start" for this number of peaks
          yy=yy-((Y2-Y1)/(X2-X1)*(xx-X1)+Y1);
          [FitResults,MeanFitError]=FitAndPlot(xx,yy,NumPeaks,Shape,delta,start,extra);
    case {49,50,51,52,53,54}
        % When a number key is pressed, sets the number of peaks to that number.  
        n=key-48;
        ShapeString=''; 
        if round(n)~=NumPeaks,
            NumPeaks=round(n);
            switch Shape
              case 1
                ShapeString='Gaussian';
              case 2
                ShapeString='Lorentzian';    
              case 3
                ShapeString='logistic';
              case 4
                ShapeString='Pearson7';
              case 5        
                ShapeString='ExpGaussian';
              otherwise
                ShapeString='';
            end % switch
          subplot(2,1,2)
          xlabel(['Number of peaks = ' num2str(NumPeaks) '    Shape = ' ShapeString  ] )
        end % if
        n=max(xx)-min(xx);
        % Create a start value for this number of peaks
        start=[];
        startpos=[n/(NumPeaks+1):n/(NumPeaks+1):n-(n/(NumPeaks+1))]+min(xx);
        for marker=1:NumPeaks,
            markx=startpos(marker);
            start=[start markx n/(5*NumPeaks)];
        end
        RedrawSignal(X,Y,xo,dx,NumPeaks);
          xlabel(['Number of peaks = ' num2str(NumPeaks) '    Shape = ' ShapeString  ] )
    case {103,108,111,112,101}
         % Selects peak shape when G,L,O,P or E kay is pressed
        switch key
            case 103 % When 'g' key is pressed, peak shape is set to Gaussian. 
                n=1;
            case 108 % When 'l' key is pressed, peak shape is set to Lorentzian. 
                n=2;
            case 111 % When 'o' key is pressed, peak shape is set to Logistic. 
                n=3;
            case 112 % When 'p' key is pressed, peak shape is set to Pearson.
                n=4;
            case 101 % When 'e' key is pressed, peak shape is set to Exponentally-broadened gaussian.
                n=5;
        end % switch
        if round(n)~=Shape,
        Shape=round(n);
          switch Shape
            case 1
                ShapeString='Gaussian';
            case 2
                ShapeString='Lorentzian';    
            case 3
                ShapeString='logistic';
            case 4
                ShapeString='Pearson';
            case 5        
                ShapeString='ExpGaussian';
            otherwise
          end % switch
          subplot(2,1,2)
          xlabel(['Number of peaks = ' num2str(NumPeaks) '    Shape = ' ShapeString  ] )
       end % if
    case 97
        % When 'a' key is pressed, increases "extra" by 0.1
        extra=extra+.1;
        if extra==0, extra=.01;,end
        [FitResults,MeanFitError]=FitAndPlot(xx,yy,NumPeaks,Shape,delta,start,extra);
    case 122
        % When 'z' key is pressed, decreases "extra" by 0.1
        extra=extra-.1;
        if extra==0, extra=.01;,end
        [FitResults,MeagfnFitError]=FitAndPlot(xx,yy,NumPeaks,Shape,delta,start,extra);
    case 114
        % When 'r' key is pressed, prints out table of fit results
       disp('Peak#  Position       Height        Width         Area') 
         for peak=1:NumPeaks,
               disp(num2str(FitResults(peak,:)))
         end % for
    case 100
        % When 'd' key is pressed, prints out table of data segment
       disp('     x          y') 
         disp(num2str([xx(1:length(xx))' yy(1:length(xx))']))
    case 116
        % When 't' key is pressed, toggles AUTOZERO mode
       AUTOZERO=-AUTOZERO    
    case 113
        % When 'q' key is pressed, prints out fitting parameters
        switch Shape
           case 1
               ShapeString='Gaussian';
           case 2
               ShapeString='Lorentzian';    
           case 3
               ShapeString='Logistic';
           case 4
               ShapeString='Pearson';
           case 5        
               ShapeString='Exponentially-broadened Gaussian';
           otherwise
         end % switch
       disp('--------------------------------------------------')
       disp(['Peak Shape = ' ShapeString])
       disp(['Number of peaks = ' num2str(NumPeaks)])
       if Shape>=4, disp(['Extra = ' num2str(extra)]), end
       disp(['Fitted range = ' num2str(min(xx)) ' - ' num2str(max(xx)) '  (' num2str(max(xx)-min(xx)) ')' ])
       disp(['Percent Error = ' num2str(MeanFitError)])
       disp('Peak#  Position       Height        Width         Area') 
       for peak=1:NumPeaks,
           disp(num2str(FitResults(peak,:)))
       end % for
   otherwise  UnassignedKey=double(key)
   end % switch
end % if
% ----------------------------------------------------------------------    
function [xx,yy,start]=RedrawSignal(X,Y,xo,dx,NumPeaks);
% Plots the entire signal (X,Y) in the lower half of the plot window and an
% isolated segment (xx,yy) in the upper half, controlled by Pan and Zoom