- 积分
- 79
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2017-3-15
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
function [smoothed,mse] = lowpass(indata,frequency,iconb,icone)(这里要怎么导入数据?)为何我的matlab运行不了?
%
% [smoothed,mse] = lowpass(indata,frequency,iconb,icone)
% lowpass filter.
% assumes unit time spaced series
%
% lowpass cutoff at f="frequency" (in cycles/time unit)
%
% icon is any one of several possible boundary constraints on the smooth at
% the begin/end x limits
% (0) minimum norm,
% (1) minimum slope
% (2) minimum roughness
%
% The lowpass routine employs a 10 point butterworth filter.
% Rather than implementing contraints (0)-(2) in the frequency domain (as in Park, 1992,
% Ghil et al, 2003) we use the following approximate implementations of the boundary
% constraints
%
% (0) pad series with long-term mean value beyond x boundary
% (1) pad series with values over last 1/2 filter width reflected w.r.t. x
% [imposes a local maximum/minimum at x boundary]
% (2) pad series with values over last 1/2 filter width reflected w.r.t. x and y (w.r.t. final value)
% [imposes a point of inflection at x boundary]
%
ipts=10;
fn=frequency*2;
% pad far enough to remove any influence of boundaries of padded series on the interior
npad=3*round(1/fn);
nn=length(indata);
padded(npad+1:npad+nn)=indata;
padded(npad+nn+1:nn+2*npad)=indata(nn);
padded(1:npad)=indata(1);
for j=nn+npad+1:nn+2*npad
ipad=j-nn-npad;
if (icone==0)
apad=mean(indata);
else if (icone==1)
apad=indata(nn-ipad);
else
apad=2*indata(nn)-indata(nn-ipad);
end
end
padded(j:j)=apad;
end
for j=1:npad
ipad=j;
if (iconb==0)
apad=mean(indata);
else if (iconb==1)
apad=indata(npad-ipad+1);
else
apad=2*indata(1)-indata(npad-ipad+1);
end
end
padded(j:j)=apad;
end
%
% smoothed=padded;
[b,a]=butter(ipts,fn,'low');
smoothed0=filtfilt(b,a,padded);
smoothed=smoothed0(npad+1:nn+npad)';
%
% determine mse relative to raw data
%
asum = 0.0;
resid=smoothed-indata;
mse=var(resid)/var(indata);
|
|