数字信号处理3个作业-----作业1带阻滤波器过滤1000HZ噪声

1. 题目

       设计一个带阻滤波器，滤除受到1000Hz正弦噪声污染的音频信号。加入正弦噪声，幅度0.1，频率1000Hz。

2.分析

      1000HZ正弦噪声的频谱图是一根(对称两根)谱线,带阻滤波器的阻带覆盖该谱线,然后变换为时域形成一个序列,将该序列与音频信号序列进行卷积即可滤去1000HZ正弦噪声,但也会滤去原信号1000HZ左右频率信号。
       使用切比雪夫1型滤波器来设计滤波器，重点是阻带通带频率的调整。然后使用切比雪夫1型滤波器获得的两个向量参数对音频信号进行卷积滤波，使用移动平均滤波器来对含噪数据进行平滑处理。此示例使用 filter 函数计算沿数据向量的平均值。移动平均值滤波器沿数据移动长度为 windowSize 的窗口，并计算每个窗口中包含的数据的平均值。以下差分方程定义向量 x 的移动平均值滤波器：

3.步骤

1.读取音频文件，格式为wav:
[audio,fs]=audioread(‘audio.wav’);
%取一个通道的信号
audio = audio(:,1); %双通道变单通道

2.产生一个正弦信号A=0.1，f=1000Hz:
noise=0.1sin(10002pi/fstt’);%加噪声 1000HZ
%tt’是(1:n)’，有公式 $\omega = 2\pi f/{f_s}$ω=2πf/fs​计算
3.选取合适的技术参数，设计滤波器，得到滤波系数[b，a]
选择带阻切比雪1型滤波器，
通带截止频率:Fp=[306 337]; %通带边界频率
通带衰减: rp=1dB
阻带截止频率:Fs=[316 327]; %阻带边界频率
阻带衰减: rs=50dB；
对应角频率(模拟滤波器的边界频率):
wp = 2*pi*Fp/Ft;
ws = 2*pi*Fs/Ft ; %求出待设计的模拟滤波器的边界频率
切比雪阶数:
[n,wn]=cheb1ord(wp,ws,rp,rs); %滤波器的最小阶数为n，wn为系统频带
滤波参数计算:
[bz,az]=cheby1(n,rp,wp,‘stop’);

4.进行滤波
z=filter(bz,az,s);

4. 结果

                                                 带阻滤波器频谱图

                                原始信号、加噪信号、带阻滤波信号比较图

5. 程序

clear
clc
[audio,fs]=audioread('audio.wav');%声音读取
audio = audio(:,1); %双通道变单通道
n=length(audio);

T = 1/fs;%采样间隔
t = (0:n-1)*T;%时间轴
f = (0:n-1)/n*fs;%频率轴

 %快速傅里叶变换
audio_fft=fft(audio,n)*T; 

%加噪声
tt =(1:n);
noise=0.1*sin(1000*2*pi/fs*tt');%加噪声 1000HZ
s=audio+noise;
s_fft=fft(s,n)*T; 
 

%设计IIR低通滤波器
rp = 1;  %通带最大衰减1db
rs=50;  %阻带最小衰减50db
Ft=fs;
Fp=[306 336]; %通带边界频率  
Fs=[316 326]; %阻带边界频率   
% Fp=4000; %通带边界频率
% Fs=5000; %阻带边界频率 
wp = 2*pi*Fp/Ft;
ws = 2*pi*Fs/Ft ;   %求出待设计的模拟滤波器的边界频率
% [N,wn]=buttord(wp,ws,rp,rs,'s');    %低通滤波器的阶数和截止频率
% [b,a]=butter(N,wn,'s');             %S域频率响应的参数即：滤波器的传输函数
% [bz,az]=bilinear(b,a,0.5);          %利用双线性变换实现频率响应S域到Z域的变换
[n,wn]=cheb1ord(wp,ws,rp,rs); %滤波器的最小阶数为n，wn为系统频带
[bz,az]=cheby1(n,rp,wp,'stop');
figure(2);%带阻滤波器特性
[h,w]=freqz(bz,az);
title('IIR带阻滤波器');
plot(w*fs/(2*pi),abs(h));  %确保1000HZ落在截断处
grid;

%滤波
z=filter(bz,az,s);
z_fft=fft(z)*T;     %滤波后的信号频谱
figure(1); 
%绘出原始音频时域波
subplot(2,3,1);
plot(t,audio);   
xlabel('时间/s');
ylabel('幅度');
title('初始信号波形');  
grid;
 %绘出原始音频频域频谱

subplot(2,3,4); 
audiof = abs(audio_fft);
plot(f(1:(n-1)/2),audiof(1:(n-1)/2));
title('初始信号频谱');
xlabel('频率/Hz');
ylabel('幅度');
grid;
%绘出加噪音频时域波
subplot(2,3,2)
plot(t,s);
title('加噪声后信号波形')
xlabel('时间/s');
ylabel('幅度');
grid;
 %绘出加噪音频频域频谱
subplot(2,3,5)
sf = abs(s_fft);
plot(f(1:(n-1)/2),sf(1:(n-1)/2));
xlabel('频率/Hz');
ylabel('幅度');
title('加噪声后信号信号频谱');
grid;
%绘出滤波音频时域波
subplot(2,3,3);
plot(t,z);
title('带阻滤波后的信号波形');
xlabel('时间/s');
ylabel('幅度');
grid;
%绘出滤波音频频域波
subplot(2,3,6);
zf = abs(z_fft);
plot(f(1:(n-1)/2),zf(1:(n-1)/2));
title('带阻滤波后信号的频谱');
xlabel('频率/Hz');
ylabel('幅度');
grid;
audio_final = [audio;s;z];%原始语音，加噪语音，滤波语音的合成音频矩阵
sound(audio_final,fs); %播放语音