%#
%# 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
自用程序:各种光谱数据预处理代码matlab.zip_EXPSMOOT_光谱数据_光谱预处理_数据预处理_预处理程序
版权申诉
5星 · 超过95%的资源 165 浏览量
2022-07-15
21:19:07
上传
评论 6
收藏 11.89MB ZIP 举报
局外狗
- 粉丝: 64
- 资源: 1万+
最新资源
- 使用 C 语言实现的计算非负整数的阶乘
- 2011-2021最新版本北京大学数字普惠金融指数(PKU-DFIIC).xlsx
- 县域数字乡村指数2018-2020(1).xlsx
- Docker容器配置进阶
- tensorflow-gpu-2.7.4-cp37-cp37m-manylinux2010-x86-64.whl
- 多段线、 圆、弧转多段线(仅我可见)
- tensorflow-2.7.2-cp38-cp38-manylinux2010-x86-64.whl
- 李慧琴C语言基础部分.zip
- yeyue-p8Yi4-ve4a83792.apk
- tensorflow-gpu-2.7.3-cp38-cp38-manylinux2010-x86-64.whl
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
评论15