数字信号处理课程设计语音信号的采样加噪滤噪还原（报告和MATLAB程序）

本设计利用计算机Windows下的录音机录入一句语音信号，然后在Matlab软件平台下，利用函数audioread对语音信号进行采样，记住采样频率和采样点数，后利用函数FFT对信号进行快速傅里叶变换，得到信号的频谱特性，然后加入一固定频率干扰信号，再画出语音信号干扰前后的时域波形，并对其频谱进行分析。最后采用双线性变换法设计几种类型的数字滤波器并对混合语音信号进行滤波、FFT快速傅里叶变换并分析各种滤波器的特点及优劣性。

clear;close all;clc;
%================原始语音信号===========================%
figure（1）;
[musicFs] = audioread（‘pei.wav‘）;                   %模数转换 对语音信号进行采样
y1=music（:1）;                                       %取单声道
y2=music（:2）;
% sound（y1Fs）;
Ns=length（y1）;                                       %采样点数Ns
t=（0:Ns-1）/Fs;                       
subplot（231）;
plot（ty1）;title（‘原始音频信号波形‘）;xlabel（‘time/s‘）;ylabel（‘幅度‘）;

Y1=fft（y1）;                                          %对y1做快速傅里叶变换
magY1=abs（Y1）;                                       %取绝对值 即模值
magY1（1）=magY1（1）/Ns;                                %0频率点的信号幅值=fft变换后的幅值/Ns
magY1（2:Ns）=magY1（2:Ns）/Ns*2;                        %其他频率点信号的幅值=fft变换后的幅值*2/Ns
subplot（234）;
f=（0:Ns-1）*Fs/Ns;                                    %Fs/Ns为频率分辨率 fft的每一点乘以分辨率代表此点的频率
plot（fmagY1）;title（‘幅频特性曲线‘）;xlabel（‘frequency/Hz‘）;ylabel（‘幅度‘）;axis（[0 1e+4  0 0.03]）;

%================噪声====================================%
t=（0:Ns-1）/Fs;             
noise=sin（2*pi*7500*t）;                              %噪声为正弦波 幅值为1 频率为7500Hz
subplot（232）;
plot（tnoise）;title（‘噪声信号波形‘）;xlabel（‘time/s‘）;ylabel（‘幅度‘）;

Noise=fft（noise）;
magN=abs（Noise）;
magN（1）=magN（1）/Ns;
magN（2:Ns）=magN（2:Ns）/Ns*2;
subplot（235）;
f=（0:length（Noise）-1）*Fs/Ns;
plot（fmagN）;title（‘幅频特性曲线‘）;xlabel（‘frequency/Hz‘）;ylabel（‘幅度‘）;axis（[0 1e+4  0 0.03]）;

%===============加性噪声================================%
y1n=y1+noise‘;                                      
% sound（y1nFs）;
t=（0:length（y1n）-1）/Fs;
subplot（233）;
plot（ty1n）;title（‘加噪后音频信号波形‘）;xlabel（‘time/s‘）;ylabel（‘幅度‘）;

Xn=fft（y1n）;
magXn=abs（Xn）;
magXn（1）=magXn（1）/Ns;
magXn（2:Ns）=magXn（2:Ns）*2/Ns;
subplot（236）;
f=（0:length（Xn）-1）*Fs/length（Xn）;
plot（fmagXn）;title（‘幅频特性曲线‘）;xlabel（‘frequency/Hz‘）;ylabel（‘幅度‘）;axis（[0 1e+4  0 0.03]）;


%==================巴特沃斯低通滤波器======================%
fp=4900;fs=7450;As=100;Rp=1;                               %缺点：过渡带要求较宽 阶数高 26阶
OmegaP=2*pi*fp;OmegaS=2*pi*fs;                             %模拟角频率
wp=OmegaP/Fs;ws=OmegaS/Fs;                                 %数字频率
wp1=2*Fs*tan（wp/2）;                                        %频率预矫正
ws1=2*Fs*tan（ws/2）;
[Nwc]=buttord（wp1ws1RpAs‘s‘）;                         %切比雪夫II型阶数与截止频率
[ba]=butter（Nwc‘s‘）;                                    %模拟切比雪夫滤波器分子、分母系数
[bzaz]=bilinear（baFs）;                                  %双线性变换AF到DF
w0=[wpws];
Hx1=freqz（bzazw0）;                                        %计算在矢量w0中指定频率处的得频率响应
[Hw]=freqz（bzaz）;
dbHx1=-20*log10（abs（Hx1）/max（abs（H）））;                       %求滤波器截止频率处的衰减
figure（4）;
plot（w/（2*pi）*Fs/100020*log10（abs（H）））;
title（‘巴特沃斯低通数字滤波器幅频特性‘）;
xlabel（‘f（KHz）‘）;ylabel（‘dB‘）;

%----滤波后音频信号----%
figure（2）;
out=filter（bzazy1n）;
t=（0:length（out）-1）/Fs;
subplot（221）;
plot（tout）;title（‘巴特沃斯低通滤波后音频信号波形‘）;xlabel（‘time/s‘）;ylabel（‘幅度‘）;

OUT=fft（out）;
magOUT=abs（OUT）;
magOUT（1）=magOUT（1）/Ns;
magOUT（2:Ns）=magOUT（2:Ns）/Ns*2;
subplot（223）;
f=（0:length（OUT）-1）*Fs/length（OUT）;
plot（fmagOUT）;title（‘幅频特性曲线‘）;xlabel（‘freq

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      61440  2016-12-21 20:22  16171课程设计资料\16171课程设计资料\数字信号处理课程设计教学日历.doc

     文件      86528  2015-12-24 16:57  16171课程设计资料\16171课程设计资料\数字信号课设教案.doc

     文件     222720  2015-12-24 17:03  16171课程设计资料\16171课程设计资料\语音信号滤波课程设计指导书.doc

     目录          0  2016-12-26 09:30  16171课程设计资料\16171课程设计资料

     目录          0  2016-12-26 09:30  16171课程设计资料

----------- ---------  ---------- -----  ----

               370688                    5