DSP单频率噪音信号处理MATLAB代码

close all;
[x,fs,N] =wavread('E:\MATLAB\Eagles - Hotel California.wav')   ;   %读声音文件
n=length(x);
t=0:1/fs:(length(x)-1)/fs;        %求出语音信号的长度
y1=fft(x,n)    ;                  %傅里叶变换  
y2=fftshift(y1);                  %对频谱图进行平移 
f=0:fs/n:fs*(n-1)/n  ;            %得出频点   
figure(1) ; 
subplot(2,1,1); 
plot(t,x)                         %做原始语音信号的时域图形 
title('原始信号时域波形图');
subplot(2,1,2);
plot(abs(y1));
axis([0,300000,0,4000]);          %确定横纵坐标的范围
title('原始信号频谱图')
%加入噪声
A=0.01;                            %噪声幅值  
d=[A*cos(2*pi*10000*t)]';          %干扰信号的构建，构建了一个余弦函数
% figure(9)
% plot(d);
% axis([0,1000,-0.03,0.03])
x1=x(:,1)+d;                       %取原始语音信号的单声部信号然后和噪声信号相加 
%sound(x1,fs,N);                   %播放加噪后的音频
y1=fft(x1,n);    
y2=fftshift(y1);                   %对频谱图进行平移
figure(2) ; 
subplot(2,1,1);                    %创建两行一列绘图区间的第1个绘图区间 
plot(t,x1);  
title('加噪后的信号时域波形图'); 
xlabel('时间'); 

ylabel('幅值'); 
subplot(2,1,2);  
plot(abs(y1)); 
title('加噪后的信号频谱'); 
axis([0,300000,0,4000]);
xlabel('频率/Hz');
ylabel('幅值/V');
%低通滤波器的设计
wp=0.256*pi;ws=0.456*pi;            %确定滤波器的通带，阻带截止频率
omegap=2*fs*tan(0.256*pi/2);
omegas=2*fs*tan(0.456*pi/2);        %双线性变换频率的预畸
rp=1;as=100;                        %通带所允许的最大衰减，阻带所达到的最小衰减
[N,omegaC]=cheb2ord(omegap,omegas,rp,as,'s');%切比雪夫2型滤波器的阶次及3dB处的截止频率
[b,a]=cheby2(N,as,omegas,'s');      %切贝雪夫2型，b,a分别为滤波器的分子多项式与分母多项式的系数向量
[bz,az]=bilinear(b,a,fs);           %进行双线性变换，实现从模拟到数字滤波器的变换
[h,w]=freqz(bz,az,fs);              %低通滤波器的频率响应
figure(5);
subplot(3,1,1);
plot(w/pi,20*log10(abs(h)));        %画出滤波器
title('低通滤波器');
y9=filter(bz,az,x1);                %通过低通滤波器进行滤波
figure(4);
subplot(2,1,1);
plot(y9);
title( '低通滤波后之后信号时域图' );  
y10 = fft (y9); 
y11 = fftshift (y10);               %对频谱图进行评平移
subplot(2,1,2);
plot(abs(y10));
axis([0,300000,0,4000]);