自用程序：各种光谱数据预处理代码matlab

%# %# function [rx,gx]=lowp(sx,win,cutoff,acq) %# %# AIM: Filtering in Fourier domain by using windowing the Fourier %# transform of the signal. %# %# PRINCIPLE: Filtering in Fourier domain corresponds to applying a low pass %# filter in the frequency domain. The denoising procedure can be %# applied to signal in which the noise component has higher %# frequencies than the analytical signal. %# It can be applied to chromatographic and spectroscopic %# data with white noise. %# H. C. Smit, Chemom. and Intell. Lab. Syst., 8, 15, 1990. %# %# INPUT: sx: (1 x n) Signal vector. (n must be a power of 2) %# win: (1 x 1) Window type (Optional) %# win=0, Hamming window (Default) %# win=1, Kaiser window %# win=2, Square window %# cutoff: (1 x 1) Cutoff frequency value (Optional) %# Represents the last value different from zero %# which is retained in the low pass filter operation. %# It can be chosen visually. %# Default : cutoff=0 (automatically calculated on the %# bases of the power spectra) %# %# OUTPUT: rx: (1 x n) vector containing the filtered signal %# gx: (1 x n) vector containing the Fourier coefficients %# %# AUTHOR: Luisa Pasti %# Copyright(c) 1997 for ChemoAC %# FABI, Vrije Universiteit Brussel %# Laarbeeklaan 103 1090 Jette %# %# VERSION: 1.1 (28/02/1998) %# %# TEST: Kris De Braekeleer %# function [rx,gx]=lowp(sx,win,cutoff,acq); % set the optional parameters if nargin == 1 win=0; cutoff=0; acq=0; end if nargin == 2 cutoff=0; acq=0; end if nargin == 3 acq=0; end if acq==0 acq=1; end % length of the signal np=length(sx); % half of the length of the vector np2=np/2; % L is the power of two that correspond to the signal length L=log2(np); % maximum frequency fm=(1/acq)*0.5; % nyquisit frequency ny=0.5; % Fourier transform of the signal fx=fft(sx); % Power spectrum of the signal px=fx.*conj(fx); % Automatically computed cutoff value % Maximum of the signal ma=max(sx); % find the index (f2) of the signal greater then half of the maximum [f1 f2]=find((sx-ma/2+ma/10)>0); % length of the signal greater then half of the maximum value nf=length(f2); % numerical derivative of the x-axis index greater then half of the maximum value df=f2(2:1:nf)-f2(1:1:nf-1); % select the largest peaks of the signal by taking % the largest interdistances between two indices [f3 f4]=find(df>=5); ndf=length(f3); % largest peak width computed as differences between % the indices related to the part of the signal larger of half of the maximum ddf=(f3(2:1:ndf)-f3(1:1:ndf-1)); % largest peak of the signal mf=max(ddf); fdf=nonzeros(ddf-1); % narrowest peak of the signal mif=min(fdf); % the mimimum frequency in Fourier domain to be keep % is related to the narrowest peaks of the signal fq=round(1/(mif+1)*np2/ny); if fq>=np/3 nk=np/3; else nk=np/3; end % Splitting of the power spectra in different ranges % on the basis of the length of the signal: np=2^L; if L>7 kk=16; else kk=8; end % step in which devide the power spectra stp=np2/(kk); % mean energy of the signal (mean value of the power spectra) % in the high frequency domain (related to noise) mt=mean((px(nk:np/2))); % Mean energy of the signal in the different ranges of frequency for i=1:kk mn(i)=mean(px(1+(i-1)*stp:i*stp)); end % line equation: Approximation of the power spectra content % by means of line. for i=1:kk-1 % slope slp(i)=(mn(i+1)-mn(i))./(i*stp-(1+(i-1)*stp)); % x-axis values for each step x(i,:)=[(1+(i-1)*stp):i*stp]; % y-axis values for each step y(i,:)=mn(i+1)-((i*stp)-x(i,:)).*slp(i); % y-values of the line yt=[yt y(i,:)]; end if cutoff==0 % find the cutoff point as the point of the line which % represent the power spectra and the mean energy of the signal % related to the noise [f5 f6]=find(yt<mt); kf=f6(1); else % input cutoff % conversion of the cutoff value from frequency to number of point cut=round(cutoff*np2/ny); kf=cut; end % change the selected cutoff value by visual selection disp(' ') an=input('## Would you like to change the selected cutoff? (y/n) ','s'); if an=='Y'; an='y'; end if an=='y', disp(' ') disp('## Use the left mouse button to pick the end point of the frequency') disp('## Press any key to continue') pause ax=[fm/(np2):fm/np2:fm]; % display the logarithm of the power spectra of the signal px % the mean energy of the signal in high frequency mt % the line approximation of the power spectra yt % the selected cutoff (*) semilogy(ax(1:np2),px(1:np2),ax(1:np2),mt*ones(1,np2),ax(1:length(yt)),yt,ax(kf),yt(kf),'*') xlabel('frequency'); ylabel('log spectrum'); title('select the last point and then press any key'); hold on clc % input cutoff by using the mouse % zeroing the variables x1 = []; y1 = []; n = 0; % new point but = 1; if but == 1 % input the new value [xi,yi,but] = ginput(1); plot(xi,yi,'go','era','back') n = n + 1; text(xi,yi,[' ' int2str(n)],'era','back'); x1 = [x1; xi]; y1 = [y1; yi]; end disp('End of data entry') pause close % new selected cutoff kf=round(x1*np2/ny); end % change the selected cutoff value by visual selection disp(' ') an=input('## Change the axis to select the cutoff? (y/n) ','s'); if an=='Y'; an='y'; end if an=='y', disp(' ') disp(['## Options 1) All the frequency axis'; ' 2) Half of the frequencies'; ' 3) Select the frequency']); an1=input(['## Input 1 or 2 or 3 ']); if an1==1 disp('## Use the left mouse button to pick the end point of the frequencies') disp('## Press any key to continue') pause ax=[fm/(np2):fm/np2:fm]; % display the logarithm of the power spectra of the signal px % the mean energy of the signal in high frequency mt % the line approximation of the power spectra yt % the selected cutoff (*) semilogy(ax(1:np2),px(1:np2),ax(1:np2),mt*ones(1,np2),ax(1:length(yt)),yt,ax(kf),yt(kf),'*') title('select the last point and then press any key'); xlabel('frequency'); ylabel('log spectrum'); hold on clc % input cutoff by using the mouse % zeroing the variables x1 = []; y1 = []; n = 0; % new point but = 1; if but == 1 % input the new value [xi,yi,but] = ginput(1); plot(xi,yi,'go','era','back') n = n + 1; text(xi,yi,[' ' int2str(n)],'era','back'); x1 = [x1; xi]; y1 = [y1; yi]; end disp('End of data entry') pause close % new selected cutoff kf=round(x1*np2/ny); elseif an1==2 disp('## Use the left mouse button to pick the end point of the frequencies') disp('## Press any key to continue') pause ax=[fm/(np2):fm/np2:fm]; % display the logarithm of the power spectra of the signal px % the mean energy of the signal in high frequency mt % the line approximation of the power spectra yt % the selected cutoff (*) semilogy(ax(1:np2/2),px(1:np2/2),ax(1:np2/2),mt*ones(1,np2/2),ax(1:length(yt)/2),yt(1:length(yt)/2),ax(kf),yt(kf),'*') xlabel('frequency'); ylabel('log spectrum'); title('select the last point and then press any key'); hold on clc % input cutoff by using the mouse % zeroing the variables x1 = []; y1 = []; n = 0; % new point but = 1; if but == 1 % input the new value [xi,yi,but] = ginput(1); plot(xi,yi,'go','era','back') n = n