声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1175|回复: 10

[综合讨论] 滤波器求助

[复制链接]
发表于 2007-2-25 23:00 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
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.

filter

filter
回复
分享到:

使用道具 举报

发表于 2007-2-26 16:05 | 显示全部楼层
原帖由 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,说明低通滤波器起了作用。但楼主要求是除尖峰以外,滤波前后的各分量都应一样,这样便不能使用低通滤波器。建议楼主使用带陷滤波器。

评分

2

查看全部评分

 楼主| 发表于 2007-2-26 16:36 | 显示全部楼层
Thanks a lot. I am a new comer to signal processing. I will try again. Thank you.
 楼主| 发表于 2007-2-26 18:32 | 显示全部楼层
终于可以打   汉字了,真 爽 阿。就是    打出      来 的字 重   叠 ,需要一个一个给它 挪 地方。
楼上解  释的真清楚。谢谢  阿。
 楼主| 发表于 2007-2-26 20:44 | 显示全部楼层
陷 波  滤波器,是不是用      这 个命令:
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 编辑 ]
 楼主| 发表于 2007-2-26 21:01 | 显示全部楼层
大家再帮下忙阿:loveliness:
发表于 2007-2-27 17:10 | 显示全部楼层
很抱歉,我没有用过idealfilter函数,在工具箱中也没有找到(用的是MATLAB6.5)。不知该函数是在MATLAB7中带有的,还是在什么工具箱中?
过去用notch都是自编的,本论坛上也有坛友给出了一个链接,楼主有兴趣可参考一下。
http://forum.vibunion.com/forum/ ... amp;extra=page%3D19
发表于 2007-2-27 18:25 | 显示全部楼层
实际上也可以用自适应的陷波器,这也是一个很好的方法。在本论坛上曾有坛友给出了程序,可参考一下:
http://forum.vibunion.com/forum/ ... amp;extra=page%3D25
(笫5部分)
 楼主| 发表于 2007-3-1 00:07 | 显示全部楼层
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);
 楼主| 发表于 2007-3-1 00:28 | 显示全部楼层
There is no idealfilter function in DSP toolbox
 楼主| 发表于 2007-3-1 01:46 | 显示全部楼层
idealfilter is a function in tstool.
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-30 15:14 , Processed in 0.074153 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表