滤波器求助

amster

I have spent 2 weeks on how to remove the cycles in my time series.
My data was sampled once half an hour, so there are 48 data point each day. I want to remove the daily cycle from the time series. The series has 3299 points. I tried a 8-order Low pass butterworth filter to remove the daily cycle, but the the result seemed strange. According to my knowledge, by means of removing the daily cycle, the resulted spectrum should be indentical as that of the original but eliminating the peak indicating the daily cycle.
See below. Could you please what is wrong with it? Pls find the attached.

fc=48/(3299*0.5)=0.03.     ％ Wrong or right to calculate the cutoff frequency?
[b,a]=butter(8, 0.03)   
y=filter(b,a,vola)      ％ ‘vola’ is the time series to be filtered.

songzy41

原帖由 amster 于 2007-2-25 23:00 发表
fc=48/(3299*0.5)=0.03.     ％ Wrong or right to calculate the cutoff frequency?
[b,a]=butter(8, 0.03)   
y=filter(b,a,vola)      ％ ‘vola’ is the time series to be filtered.

1，一个概念问题：在计算滤波器系数中，Wn=fc/（fs/2），其中fs是采样频率（而不应该是数据长度的一半），fc是截止频率。在楼主的问题中，每半小时采一个数据，如果以Hz为单位，fs=1/1800Hz=0.000555Hz，24小时周期对应的频率为f=1/(24*3600)=0.00001157Hz。简化起见，可用（小时）^-1为单位，fs=2，f=1/24=0.04166。截止频率fc一定要低于f，例如可取fc=0.025（虽然与楼主的数值0.03很接近，但概念不同，计算方法不同）。这样可计算滤波器系数：
Wn=fc/（fs/2）=0.025;
[b,a]=butter(8,Wn);
2，从楼主提供的二张图上看，在笫1张图中，归一化频率0.03处有一个尖峰，而笫2张图中0.03处比笫张图衰减了25dB，说明低通滤波器起了作用。但楼主要求是除尖峰以外，滤波前后的各分量都应一样，这样便不能使用低通滤波器。建议楼主使用带陷滤波器。

amster

Thanks a lot. I am a new comer to signal processing. I will try again. Thank you.

amster

终于可以打   汉字了，真 爽 阿。就是    打出      来 的字 重   叠 ，需要一个一个给它 挪 地方。
楼上解  释的真清楚。谢谢  阿。

amster

陷 波  滤波器，是不是用      这 个命令：
y= idealfilter(x,interval,'notch');
其中这个interval 怎么确定？
我这个案例里，cutoff frequency  fc=0.04166, 是不是这个interval只要把  0.04166包在里面就可以了？     如interval=[0.040,0.045]?

还是用digital signal processing 里面的bandstop？

多解指教

[ 本帖最后由 amster 于 2007-2-26 20:59 编辑 ]

amster

大家再帮下忙阿:loveliness:

songzy41

很抱歉，我没有用过idealfilter函数，在工具箱中也没有找到（用的是MATLAB6.5）。不知该函数是在MATLAB7中带有的，还是在什么工具箱中？
过去用notch都是自编的，本论坛上也有坛友给出了一个链接，楼主有兴趣可参考一下。
http://forum.vibunion.com/forum/ ... amp;extra=page%3D19

songzy41

实际上也可以用自适应的陷波器，这也是一个很好的方法。在本论坛上曾有坛友给出了程序，可参考一下：
http://forum.vibunion.com/forum/ ... amp;extra=page%3D25
（笫5部分）

amster

N=400; %总采样长度
t=0:N-1; %时间的变化范围
s=sin(2*pi*t/20); %输入信号
A=0.5; %干扰信号的幅值
fai=pi/3;%干扰信号的相移
n=A*cos(2*pi*t/10+fai);%干扰信号
x=s+n;%信号混合
subplot(2,2,1);%作第一子图
plot(t,s);
subplot(2,2,2); %作第二子图
plot(t,x);
x1=cos(2*pi*t/10);
x2=sin(2*pi*t/10);
%初始化
w1=0.1;
w2=0.1;
e=zeros(1,N);
y=0;
u=0.05;%迭代步长
for i=1:N
y=w1*x1(i)+w2*x2(i);
e(i)=x(i)-y;%误差信号
w1=w1+u*e(i)*x1(i);%迭代方程
w2=w2+u*e(i)*x2(i);%迭代方程
end
subplot(2,2,3); %作第三子图
plot(t,e);
subplot(2,2,4); %作第四子图
plot(t,s-e);

amster

There is no idealfilter function in DSP toolbox

amster

idealfilter is a function in tstool.