基于FIR的音乐加噪降噪及频谱分析

具体关于FIR滤波的原理相关的内容可以看一下其他都有详细的介绍，这里仅给出程序有详细的备注信息，能够看懂。

 

 

% 程序名称：c.m
% 程序功能：利用FLATTOPWIN设计的FIR滤波对语音信号进行滤波去噪
%audioread()
[x,fs]=audioread('Lemon.wav');  % 输了参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。
sound(x,fs);      % 按指定的采样率和每样本编码位数回放
N=length(x);       % 计算信号X的长度
fn=2200;            % 单频噪声频率
t=0:1/fs:(N-1)/fs;         % 计算时间范围，样本数除以采样频率
x=x(:,1)';               % 将双声道转为单声道
y=x+0.1*sin(fn*2*pi*t);     % 加噪声
sound(y,fs);         % 应该可以明显听出有尖锐的单品啸叫声
X=abs(fft(x));  Y=abs(fft(y));    % 对原始信号和加噪声信号进行fft变换，取幅度值
X=X(1:N/2); Y=Y(1:N/2);       % 截取前半部分
deltaf=fs/N;            % 计算频谱的谱线间隔
f=0:deltaf:fs/2-deltaf;     % 计算频谱频率范围
figure(1)

subplot(2,2,1);plot(t,x);xlabel('时间t');
ylabel('幅度');title('原始语音信号');
axis([0,4,-1.5,1.5]);

subplot(2,2,2);plot(f,X);xlabel('频率f');ylabel('幅度谱');
title('原始语音信号幅度谱');
axis([-10,6000,0,700]);
subplot(2,2,3);plot(t,y);xlabel('时间');ylabel('幅度');
title('加干扰后的语音信号');
axis([0,4,-1.5,1.5]);
subplot(2,2,4);plot(f,Y);xlabel('频率 f');ylabel('幅度谱');
title('加干扰后的语音信号幅度谱');
axis([-10,6000,0,700]);
 
fpd=2100;fsd=2150;fsu=2250;fpu=2300;Rp=1;As=30;      % 阻带滤波器设计指标
fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min((fsd-fpd),(fpu-fsu));  % 计算上下边带中心频率和频率间隔
wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;    %将Hz为单位的模拟频率换算为rad为单位的数字频率 
wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;
M=ceil(6.2*pi/dw)+1;      %计算flattopwin窗设计滤波器时需要的阶数
n=0:M-1;         % 定义时间范围
w_fla=(flattopwin(M));    % 产生M阶的flattopwin窗
hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); 
% 调用自编函数计算理想阻带滤波器的脉冲响应
h_bs=w_fla'.*hd_bs;       % 用窗口法计算实际滤波器脉冲响应
[db,mag,pha,grd,w]=freqz_m(h_bs,1);  % 调用自编函数计算滤波器的频率特征
 
figure(2)
subplot(2,2,1);plot(w,db);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('db');
axis([0,0.8,-60,10]);
line([0,2],[-As,-As],'color','y','linestyle','--','LineWidth',2);
line([0,2],[-Rp,-Rp],'color','g','linestyle','--','LineWidth',2);
line([wsd,wsd],[-30,10],'color','r','linestyle','--','LineWidth',2);
line([wsu,wsu],[-30,10],'color','r','linestyle','--','LineWidth',2);
subplot(2,2,2);plot(w,mag);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('幅度mag');
axis([0,1,-0.5,1.5]);
subplot(2,2,3);plot(w,pha);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('相位pha');
axis([0,1,-4,4]);
subplot(2,2,4);stem(n,h_bs);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('db');
axis([0,1500,0,1]);
 
y_fil=fftfilt(h_bs,y);  %用设计好的滤波器对y进行滤波
Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2);  %计算频谱取前一半
figure(3)
subplot(3,2,1);plot(t,x);xlabel('时间t');ylabel('幅度');


title('原始语音信号');
subplot(3,2,2);plot(f,X);xlabel('频率f');ylabel('幅度谱');
title('原始语音信号幅度谱');
axis([0,4000,0,600]);
subplot(3,2,3);plot(t,y);xlabel('时间t');ylabel('·幅度');
title('加干扰后语音信号');
subplot(3,2,4);plot(f,Y);xlabel('频率f');ylabel('幅度谱');
title('加干扰后语音信号幅度谱');
axis([0,4000,0,600]);
subplot(3,2,5);plot(t,y_fil);xlabel('时间t');ylabel('幅度');
title('滤波后语音信号');
subplot(3,2,6);plot(f,Y_fil);xlabel('频率f');ylabel('幅度谱');
title('滤波后语音信号幅度谱');
axis([0,4000,0,600]);
figure(4)
subplot(2,1,1);plot(f,20*log10(Y));grid on; 
subplot(2,1,2);plot(f,20*log10(Y_fil));grid on; 
sound(x,fs);   % 播放原始信号
%audiowrite('Lemon_noise.wav',x,fs);
sound (y_fil,fs);  %播放滤波后的信号
 
[C,B,A]=dir2fs(h_bs);

实验结果