查看: 899|回复: 7|关注: 0

[已解决] matlab实现滤波

[复制链接]

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
                              
有正弦信号,包含了两个频率的信息,设计一个滤波器,滤掉其中一个频率,比如x=sin(2*pi*f1*t)+sin(2*pi*f2*t);
滤掉f2的部分,我设计了一个低通滤波,f1可以通过,对x进行fft变换,得到X,X和低通滤波相乘再iff,得到f1的波形 但是幅值减半了,很不解。

怎么实现,希望大牛给与指点?
附上我的代码:
fs=1000;
N = 1024;    % fft点数
T=1/fs;w0=20;
w=0.2;
t = -w:T:w;
y1 = cos(2*pi*w0*t);
y2=cos(2*pi*5*w0*t);
y=y1+y2;
H_w=[ones(1, 2*w0) zeros(1,N-2*w0)]; %the filer freqency
% 原始信号
figure(1);
subplot(211);plot(t,y1);title('the first component');
subplot(212);plot(t,y);grid on;title('y');
%fft后的频谱图
yfft=fft(y,N);        %对yN点fft
y_w =H_w.*yfft;    %加入滤波器,时域卷积等于频域相乘;
yfftabs=abs(yfft);  
y_wabs=abs(y_w);
figure(2);
wfft=(1:N)/N*fs;
subplot(211);
plot(wfft(1:N/2),yfftabs(1:N/2));
subplot(212);
plot(wfft(1:N/2),y_wabs(1:N/2),'r');
%%傅里叶反变换
yfft_in=ifft(yfft);
yw_in=ifft(y_w);
figure(3);
subplot(211);
plot(yfft_in);grid on;
subplot(212);
plot(real(yw_in),'r');

入门

55 麦片

财富积分


50500


0

主题

71

帖子

9

最佳答案
发表于 2018-11-20 11:25:48 | 显示全部楼层
对称的,让你吃了一半

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2018-11-20 11:34:11 | 显示全部楼层
zezedi 发表于 2018-11-20 11:25
对称的,让你吃了一半

您好,plot(real(yw_in),'r');最后一句为啥要加real,这个对称是什么时候出现的,对yfft进行iff变换幅值并没有影响,能够画出原始信号,希望您能够给与解答,感谢牛人

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2018-11-20 11:45:14 | 显示全部楼层
zezedi 发表于 2018-11-20 11:25
对称的,让你吃了一半

我应该明白了你说的吃了一半,原始信号的频谱是对称的,所以用我写的这个Hw去相乘只保留了一半的频谱,是这个意思吧,我看网上都说用什么巴沃特滤波器,我这个也可以设计成那样的是吧?

入门

55 麦片

财富积分


50500


0

主题

71

帖子

9

最佳答案
发表于 2018-11-20 11:52:39 | 显示全部楼层 |此回复为最佳答案
lijie19850112 发表于 2018-11-20 11:45
我应该明白了你说的吃了一半,原始信号的频谱是对称的,所以用我写的这个Hw去相乘只保留了一半的频谱,是 ...

嗯,你这种直接在频域截断再逆变换的方法,在工程上一般很难实现,正常都是在时域卷积一个序列,这个序列就是要设计的滤波器。滤波器的频响就是你的设计要求

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2018-11-20 13:14:22 | 显示全部楼层
zezedi 发表于 2018-11-20 11:52
嗯,你这种直接在频域截断再逆变换的方法,在工程上一般很难实现,正常都是在时域卷积一个序列,这个序列 ...

谢谢,我去研究研究设计成低通滤波器的方法

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2018-11-21 20:18:59 | 显示全部楼层
lijie19850112 发表于 2018-11-20 13:14
谢谢,我去研究研究设计成低通滤波器的方法

clc;
close all;
clear all;
fs=1000;
N = 1024;    % fft点数
T=1/fs;fl=20;fh=100;
w=0.2;
t = -w:T:w;
y1 = cos(2*pi*fl*t);
y2=cos(2*pi*fh*t);
y=y1+y2;     %含有两个频率的信号
% 原始信号
figure(1);
subplot(211);plot(t,y1);title('the first component');
subplot(212);plot(t,y);grid on;title('y');
%fft后的频谱图
yfft=fft(y,N);        %对yN点fft
yfftabs=abs(yfft);  
figure(2);
wfft=(1:N)/N*fs;
plot(wfft(1:N/2),yfftabs(1:N/2));
%设计巴特沃斯滤波器
wp = 40/fs/2;    %通带截止频率
ws= 800/fs/2;      % 阻带截止频率   这两行都要处以2是为什么?
[n,Wn]=buttord(wp,ws,3,40);
[a,b]=butter(n,Wn);
    [h,f]=freqz(a,b,'whole',fs);        %求数字低通滤波器的频率响应
    f=(0:length(f)-1*fs/length(f));     %进行对应的频率转换
    figure;
    subplot(311);
    plot(f(1:length(f)/2),abs(h(1:length(f)/2)));       %m 绘制幅频响应图
    title('Butterworth low-pass filter ');
    xlabel('frequency/Hz');
    ylabel('Amplitude');grid;
   
    %对有两个频率的信号进行滤波
    sF=filter(a,b,y);  
%     figure;
    subplot(312);
    plot(sF);                      %plot signal after filter
    title('滤波之后的时域信号 ws=800/fs/2');xlabel('t/s');ylabel('Amplitude');grid;
    SF=fft(sF,1024);  %滤波后的信号进行fft,看看你频谱中包含的频谱分量
    subplot(313);
%     plot((1:length(SF)/2)*fs/length(SF),2*abs(SF(1:length(SF)/2))/length(SF));
plot(wfft(1:N/3),2*abs(SF(1:N/3)));
    title('滤波之后的信号的频谱');grid;
    xlabel('frequency/Hz');
    ylabel('Amplitude');
您好,这是我采用的巴特沃斯滤波器,wp和ws到底要怎么设置,仿真的图形差别很大,而且有失真,望牛人给与指点

ws1.png
ws2.png

新手

5 麦片

财富积分


050


1

主题

16

帖子

0

最佳答案
发表于 2018-12-13 19:02:50 | 显示全部楼层
你好,请问最后解决了吗?能分享一下吗?多谢啦
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

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