FIR 滤波器-凯泽(Kaiser)窗 matlab源程序

clc; % 清除工作区 clear; % 清除命令窗口 close all; % 关闭窗口 %% 读取音频文件 [x,fs]=audioread('audio.mp3'); % 读取音频信号：x是数据，fs是采样率 x = x(:,1)'; % 取第一列数据转置 N=length(x); % 计算信号x的长度 t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率 % sound(x,fs); % 播放音频 X=abs(fft(x));X=X(1:N/2); % 对原始信号进行ff变换，截取幅度谱前半部分 deltaf=fs/N; % 计算频谱的谱线间隔 f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围 figure(1) subplot(2,1,1);plot(t,x);axis([0 6 -1 1]);grid on; % 原始语音信号的时域图 xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16); title('原始语音信号','FontSize',20); subplot(2,1,2);plot(f,X);axis([0 15000 0 1200]);grid on; % 画原始语音信号幅度谱图 xlabel('频率/Hz','FontSize',16);ylabel('|X(f)|','FontSize',16); title('原始语音信号幅度谱图','FontSize',20); %% 生成窄带高斯白噪声 noise = wgn(1,N,-20); [n,Wn,beta,ftype] = kaiserord([10000 10100 10900 11000],... [0 1 0],[0.01 0.02 0.01],fs); % 估计kaiser窗的参数 n=n+rem(n,2); % 保证滤波器系数长N+1为奇数 b_noise=fir1(n,Wn,ftype,kaiser(n+1,beta)); % kaiser窗函数滤波器设计 noise_fil=fftfilt(b_noise,noise); % 滤出窄带噪声 %% 对原始信号加窄带高斯白噪声 y=x+noise_fil; % 将窄带高斯白噪声叠加到原始音频上 % sound(y,fs); % 播放加噪信号 Y=abs(fft(y));Y=Y(1:N/2); % 对加噪信号进行ff变换，截取幅度谱前半部分 figure(2) subplot(2,1,1);plot(t,y);axis([0 6 -1 1]);grid on; % 加噪语音信号的时域图 xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16); title('加入干扰后的语音信号','FontSize',20); subplot(2,1,2);plot(f,Y);axis([0 15000 0 1200]);grid on; % 加噪语音信号频谱 xlabel('频率/Hz','FontSize',16);ylabel('|Y(f)|','FontSize',16); title('加入干扰后的语音信号幅度谱图','FontSize',20); %% 滤波器设计 Fs2=fs/2; % 奈奎斯特频率 fp=9900; % 通带频率 fr=10000; % 阻带频率 Rp=3; % 通带最大波纹衰减 As=40; % 阻带最小衰减 delta1 = (10^(Rp/20)-1)/(10^(Rp/20)+1); % 通带波纹线性值 delta2 = (1+delta1)*(10^(-As/20)); % 阻带衰减线性值 ff=[fp fr]/Fs2; % 频率指标（归一化频率） A=[1 0]; % 和幅值指标 dev=[delta1 delta2]; % 偏离指标（偏离程度） [M,Wn,beta,ftype] = kaiserord(ff,A,dev); % 估计kaiser窗的参数 M=M+rem(M,2); % 保证滤波器系数长N+1为奇数 n=0:M; % 时间范围 b=fir1(M,Wn,ftype,kaiser(M+1,beta)); % kaiser窗函数滤波器设计 [db,mag,pha,grd,w]=freqz_m(b,1); % 计算频率响应、相位响应、脉冲响应 figure(3) subplot(2,2,1);plot(w/pi,db);axis([0 1 -150 10]);grid on; % 幅度响应(DB) xlabel('w/pi','FontSize',16);ylabel('dB','FontSize',16); title('滤波器幅度响应图','FontSize',20); subplot(2,2,2);plot(w/pi,mag);axis([0 1 -0.2 1.2]);grid on; % 幅度响应 xlabel('wp','FontSize',16);ylabel('幅度 mag','FontSize',16); title('滤波器幅度响应图','FontSize',20); subplot(2,2,3);plot(w/pi,pha);grid on; % 相位响应 xlabel('w/pi','FontSize',16);ylabel('相位 pha','FontSize',16); title('滤波器相位响应图','FontSize',20); subplot(2,2,4);stem(n,b);grid on; % 脉冲响应 xlabel('n','FontSize',16);ylabel('h(n)','FontSize',16); title('滤波器脉冲响应图','FontSize',20); %% 信号滤波处理 y_fil=fftfilt(b,y); % 用设计好的滤波器对y进行滤波 Y_fil=abs(fft(y_fil)); Y_fil=Y_fil(1:N/2); % 计算频谱取前一半 figure(4) subplot(3,2,1);plot(t,x);axis([0 6 -1 1]);grid on; xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16); title('原始语音信号','FontSize',20); subplot(3,2,2);plot(f,X);axis([0 15000 0 1200]);grid on; xlabel('频率/Hz','FontSize',16);ylabel('|X|(f)','FontSize',16); title('语音信号幅度谱','FontSize',20); subplot(3,2,3);plot(t,y);axis([0 6 -1 1]);grid on; xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16); title('加干扰后的语音信号','FontSize',20); subplot(3,2,4);plot(f,Y);axis([0 15000 0 1200]);grid on; xlabel('频率/Hz','FontSize',16);ylabel('|Y|(f)','FontSize',16); title('加干扰语音信号幅度谱','FontSize',20); subplot(3,2,5);plot(t,y_fil);axis([0 6 -1 1]);grid on; xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16); title('滤波后语音信号','FontSize',20); subplot(3,2,6);plot(f,Y_fil);axis([0 15000 0 1200]);grid on; xlabel('频率/Hz','FontSize',16);ylabel('|Y\_fil|(f)','FontSize',16); title('滤波后语音信号幅度谱','FontSize',20); % sound(y_fil,fs); % 播放音频