matlab给音频信号加噪并设计IIR滤波器滤去噪声

先记录一些MATLAB的常用指令，方便看后面代码时参考。

matlab的一些常用指令

matlab官方网站上给出了所有指令的详细文档，这里只列出一丢丢代码里用到的：），需要用到时去上面查就可以啦。

• audioread()

[y,Fs] = audioread(filename) 从名为 filename 的文件中读取数据，并返回样本数据 y 以及该数据的采样率 Fs。

• fft()

Y = fft(X,n) 返回 n 点 DFT。如果未指定任何值，则 Y 的大小与 X 相同。

• 如果 X 是向量且 X 的长度小于 n，则为 X 补上尾零以达到长度 n。
• 如果 X 是向量且 X 的长度大于 n，则对 X 进行截断以达到长度 n。
• 如果 X 是矩阵，则每列的处理与在向量情况下相同。
• 如果 X 为多维数组，则大小不等于 1 的第一个数组维度的处理与在向量情况下相同。

• figure()

figure(n) 查找 Number 属性等于 n 的图窗，并将其作为当前图窗。如果不存在具有该属性值的图窗，MATLAB® 将创建一个新图窗并将其 Number 属性设置为 n。

• subplot()

subplot(m,n,p) 将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。MATLAB® 按行号对子图位置进行编号。第一个子图是第一行的第一列，第二个子图是第一行的第二列，依此类推。如果指定的位置已存在坐标区，则此命令会将该坐标区设为当前坐标区。

• plot()

plot(X,Y) 创建 Y 中数据对 X 中对应值的二维线图。

• 如果 X 和 Y 都是向量，则它们的长度必须相同。plot 函数绘制 Y 对 X 的图。
• 如果 X 和 Y 均为矩阵，则它们的大小必须相同。plot 函数绘制 Y 的列对 X 的列的图。
• 如果 X 或 Y 中的一个是向量而另一个是矩阵，则矩阵的各维中必须有一维与向量的长度相等。如果矩阵的行数等于向量长度，则 plot 函数绘制矩阵中的每一列对向量的图。如果矩阵的列数等于向量长度，则该函数绘制矩阵中的每一行对向量的图。如果矩阵为方阵，则该函数绘制每一列对向量的图。
• 如果 X 或 Y 之一为标量，而另一个为标量或向量，则 plot 函数会绘制离散点。但是，要查看这些点，您必须指定标记符号，例如 plot(X,Y,‘o’)。

• title()

title(txt) 将指定的标题添加到 gca 命令返回的坐标区或图中。重新发出 title 命令可使新标题替换旧标题。

• xlabel() & ylabel()

xlabel(txt) 将为 gca 命令返回的当前坐标区或图的 x 轴添加标签。重新发出 xlabel 命令会将旧标签替换为新标签。
ylabel(txt) 将为 gca 命令返回的当前坐标区或图形的 y 轴添加标签。重新发出 ylabel 命令可使新标签替换旧标签。

• sound()

sound(y,Fs) 以采样率 Fs 向扬声器发送音频信号 y。

• buttord()

• [n,Wn] = buttord(Wp,Ws,Rp,Rs)返回数字巴特沃斯滤波器的最低阶n，其通带纹波不超过Rp dB，阻带衰减至少为Rs dB。Wp和Ws分别滤波器的通带和阻带边缘频率,规范化从0到1,1对应于πrad /样品。相应截止频率的标量(或向量)Wn也被返回。要设计一个Butterworth过滤器，使用输出参数n和Wn作为butter的输入。
• [n,Wn] = buttord(Wp,Ws,Rp,Rs，‘s’)求出模拟巴特沃思滤波器的最小阶n和截止频率。指定频率Wp和Ws(以弧度/秒为单位)。通带和阻带可以是无限的。

• butter()

[b,a] = butter(n,Wn)返回具有归一化截止频率Wn的n阶低通数字巴特沃思滤波器的传递函数系数。

• bilinear()

[numd,dend] =bilinear(num,den,fs)将分子num和分母den指定的s域传递函数转换为离散的等价函数。

• freqz()

数字滤波器的频率响应
[h,w] = freqz(b,a,n)返回存储在b和a中的传递函数系数的数字滤波器的n点频响向量h和对应的角频率向量w。

• randn()

X = randn(sz1,…,szN) 返回由随机数组成的 sz1×…×szN 数组，其中 sz1,…,szN 指示每个维度的大小。例如：randn(3,4) 返回一个 3×4 的矩阵。

直接给出代码（最近考试略多，以后有时间再详细解释）

%% read audio from file
filename = 'D:\matlab_exercise\audio_proc\testaudio.wav';
[y,fs] = audioread(filename); % read audio from file, return sample data y and sampling rate fs 
%sound(y,fs); % play the audio

%% show the time domain waveform of the original signal
N = length(y); % number of samples 
time = 1/fs * (0:N-1); % 1/fs is sample interval
figure(1)
subplot(4,2,1);
plot(time,y);
title('The sampled time domain waveform of original signal');
grid on;
xlabel('Time/s');
ylabel('Amplitude');

%% show the frequency domain waveform of the original signal
y_FFT = fft(y,N); % returns the n-point DFT
f_p = fs / N * (0:N/2-1);
subplot(4,2,2);
plot(f_p,y_FFT(1:N/2));
title('The sampled frequency domain waveform of original signal');
grid on;
xlabel('f/Hz');
ylabel('X(f)');

%% add noise
noise = 0.01 * randn(N,1); % generate random noise of equal length
y_mix = y + noise; % mix original audio and noise
%sound(y_mix,fs);

%% show the time domain waveform of the mixed signal
subplot(4,2,3);
plot(time,y_mix);
title('The sampled time domain waveform of mixed signal');
grid on;
xlabel('Time/s');
ylabel('Amplitude');

%% show the frequency domain waveform of the mixed signal
y_mix_FFT = fft(y_mix,N);
subplot(4,2,4);
plot(f_p,y_mix_FFT(1:N/2));
title('The sampled frequency domain waveform of mixed signal');
grid on;
xlabel('f/Hz');
ylabel('X(f)');

%% digital low-pass technical index
Ft = fs / N;
Fp = 3000 * Ft; % 0.8*10^(4)
Fs = 3600 * Ft; % 0.85*10^(4)
Wp = 2 * pi * Fp;
Ws = 2 * pi * Fs;
wp = Wp / fs;
ws = Ws / fs;
alpha_p = 1;
alpha_s = 10;

%% analog low-pass technical index
T = 1;
omega_p = 2/T * tan(wp / 2);
omega_s = 2/T * tan(ws / 2);

%% design analog filter
[n_order,wc] = buttord(omega_p,omega_s,alpha_p,alpha_s,'s'); % wc : 3dB cut-off frequency
[B,A] = butter(n_order,wc,'s');

%% convert to digital filter
[Bz,Az] = bilinear(B,A,1/T);
[Hz,wk] = freqz(Bz,Az,N);

%% show filter
subplot(4,2,5);
plot(wk,abs(Hz));
%plot(wk/pi,20*log10(abs(Hz)));
grid on;
xlabel('\omega/\pi');
ylabel('Amplitude(dB)');
%axis([0,1,-100,5]);
title('filter');

%% filter out the noise
y_filt = filter(Bz,Az,y_mix);
sound(y_filt,fs);

%% show the time domain waveform after filter
subplot(4,2,7);
plot(time,y_filt);
title('The sampled time domain waveform of mixed signal');
grid on;
xlabel('Time/s');
ylabel('Amplitude');

%% show the frequency domain waveform after filter
y_filt_FFT = fft(y_filt,N);
subplot(4,2,8);
plot(f_p,y_filt_FFT(1:N/2));
title('The sampled frequency domain waveform of mixed signal');
grid on;
xlabel('f/Hz');
ylabel('X(f)');