skaly2008
2009-03-01, 14:38
matlab过滤器的缺点?
最近找了一个过滤器,发现滤波器还是有问题啊。。。别的不说,请看代码:
clear;
fc1=22999;%通带下限截止频率fc1=22999 Hz
fc2=23001;%通带上限截止频率fc2=23001 Hz
fs=102400;%采样频率fs=102400Hz
N= 50000;%采样点数目
t = (0:N)/fs;
f1=5000;
f2=15000;
f3=30000;
f4=23000;
f5=50000;
f6=60000;
s1=CreateLine(f1,t);
s2=CreateLine(f2,t);
s3=CreateLine(f3,t);
s4=CreateLine(f4,t);
s5=CreateLine(f5,t);
s6=CreateLine(f6,t);
s=s1+s2+s3+s4+s5+s6;
sf = beltfilter(s,fs,20,fc1,fc2);%对信号s进行滤波
figure(1);
subplot(311);
plot(t,s);%作正弦信号的时域波形
xlabel('t');
ylabel('y');
title('原混合波形');
figure(1);
subplot(312);
xl=length(sf);
t=t(1:xl);
plot(t,sf);%作正弦信号的时域波形
xlabel('t');
ylabel('y');
title('过滤后波形');
subplot(313);
s3=s3(1:xl);
plot(t,s3);%作正弦信号的时域波形
xlabel('t');
ylabel('s3');
title('原波形');
N=xl;
xf = 2*abs(fftn(sf))/N;
xf = xf(1:N/2);
df = fs/N; %频率分辨率Hz
f = (0:N/2-1)*df; %频域序列
cMax = max(xf) ;
figure(2);
subplot(1,1,1);
plot(f, xf);
xlabel('f/Hz');
f(find(xf>=max(xf)))
CreateLine:
%创建一个序列,随机把一些连续值清空
%a-频率
%t-采样时间点序列
function xl=CreateLine(a,t)
N=length(t);
%振幅
iZF=fix(100*rand);
%初始相位
iCSXW=rand;
%开始清除点
iStart=fix(N*rand);
%结束清除点
iStop=fix(N*rand);
%使iStart<iStop
if iStart>iStop
iTmp=iStop;iStop=iStart;iStart=iTmp;
end;
xl=iZF*sin(a*2*pi*t+iCSXW);
for k=iStart:iStop
xl(k)=0;
end;
beltfilter:(注:网上找到的)
%带通滤波
%x-输入信号
%fs-输入信号采样频率
%M-滤波器半阶数
%fl-下截止频率
%fh-上截止频率
%xl-返回的滤波结果,其长度为原序列长度减少2M,去掉了前M和后M个不正确的点
function xl=beltfilter(x,fs,M,fl,fh)
N=length(x)-2*M;
[i j]=size(x);
if i~=1
x=x';
end
k=1:M;
w=0.5+0.5*cos(pi*k/M);
wl=2*pi*fl/fs;
wh=2*pi*fh/fs;
h(1)=(wh-wl)/pi;
h(2:M+1)=(sin(wh*k)-sin(wl*k))./(pi*k).*w;
for k=1:N
kk=k-1+M;
xl(k)=x(kk+1)*h(1)+sum(h(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
end
if i~=1
xl=xl';
end
大家看看结果就知道了$$$$
http://www.ilovematlab.cn/attachments/month_0903/20090301_9ef74156e522bb714b292eYSeP6fOhNL.jpg
最近找了一个过滤器,发现滤波器还是有问题啊。。。别的不说,请看代码:
clear;
fc1=22999;%通带下限截止频率fc1=22999 Hz
fc2=23001;%通带上限截止频率fc2=23001 Hz
fs=102400;%采样频率fs=102400Hz
N= 50000;%采样点数目
t = (0:N)/fs;
f1=5000;
f2=15000;
f3=30000;
f4=23000;
f5=50000;
f6=60000;
s1=CreateLine(f1,t);
s2=CreateLine(f2,t);
s3=CreateLine(f3,t);
s4=CreateLine(f4,t);
s5=CreateLine(f5,t);
s6=CreateLine(f6,t);
s=s1+s2+s3+s4+s5+s6;
sf = beltfilter(s,fs,20,fc1,fc2);%对信号s进行滤波
figure(1);
subplot(311);
plot(t,s);%作正弦信号的时域波形
xlabel('t');
ylabel('y');
title('原混合波形');
figure(1);
subplot(312);
xl=length(sf);
t=t(1:xl);
plot(t,sf);%作正弦信号的时域波形
xlabel('t');
ylabel('y');
title('过滤后波形');
subplot(313);
s3=s3(1:xl);
plot(t,s3);%作正弦信号的时域波形
xlabel('t');
ylabel('s3');
title('原波形');
N=xl;
xf = 2*abs(fftn(sf))/N;
xf = xf(1:N/2);
df = fs/N; %频率分辨率Hz
f = (0:N/2-1)*df; %频域序列
cMax = max(xf) ;
figure(2);
subplot(1,1,1);
plot(f, xf);
xlabel('f/Hz');
f(find(xf>=max(xf)))
CreateLine:
%创建一个序列,随机把一些连续值清空
%a-频率
%t-采样时间点序列
function xl=CreateLine(a,t)
N=length(t);
%振幅
iZF=fix(100*rand);
%初始相位
iCSXW=rand;
%开始清除点
iStart=fix(N*rand);
%结束清除点
iStop=fix(N*rand);
%使iStart<iStop
if iStart>iStop
iTmp=iStop;iStop=iStart;iStart=iTmp;
end;
xl=iZF*sin(a*2*pi*t+iCSXW);
for k=iStart:iStop
xl(k)=0;
end;
beltfilter:(注:网上找到的)
%带通滤波
%x-输入信号
%fs-输入信号采样频率
%M-滤波器半阶数
%fl-下截止频率
%fh-上截止频率
%xl-返回的滤波结果,其长度为原序列长度减少2M,去掉了前M和后M个不正确的点
function xl=beltfilter(x,fs,M,fl,fh)
N=length(x)-2*M;
[i j]=size(x);
if i~=1
x=x';
end
k=1:M;
w=0.5+0.5*cos(pi*k/M);
wl=2*pi*fl/fs;
wh=2*pi*fh/fs;
h(1)=(wh-wl)/pi;
h(2:M+1)=(sin(wh*k)-sin(wl*k))./(pi*k).*w;
for k=1:N
kk=k-1+M;
xl(k)=x(kk+1)*h(1)+sum(h(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
end
if i~=1
xl=xl';
end
大家看看结果就知道了$$$$
http://www.ilovematlab.cn/attachments/month_0903/20090301_9ef74156e522bb714b292eYSeP6fOhNL.jpg