clc;
clear all;
w=1:20000;
[y,Fs,bits] = wavread('f41.wav');
L_Time_sound=length(y)/Fs;
Time_scale=0:L_Time_sound/length(y):L_Time_sound-L_Time_sound/length(y);
Max_analog_F=Fs/2;%采样率为 Fs 的数字信号至多能表示 Fs/2 的模拟信号
Frequency_scale=-Max_analog_F:2*Max_analog_F/length(y):Max_analog_F-Max_analog_F/len
gth(y);
figure(1);
subplot(2,3,1);
plot(Time_scale,y);
title('时域波形');xlabel('时间长度/秒')
Y=fft(y);
amp_Y=fftshift(abs(Y));
subplot(2,3,4);
plot(Frequency_scale,amp_Y);
xlabel('模拟频率/Hz');ylabel('|Y|');title('频域幅度谱');
%%
% fb=300HZ,阻带频率 fc=1200HZ,通带波纹 Ap=1dB,阻带波纹 As=40dB
%
fp=500;ffs=1000;Rp=1;Rs=40;T=1/Max_analog_F;
W1p=fp/Max_analog_F*2; W1s=ffs/Max_analog_F*2;
[N, Wn] = buttord(W1p, W1s, Rp, Rs, 's');
[z,p,k] = buttap(N);
[bp,ap] = zp2tf(z,p,k);
[bs,as] = lp2lp(bp,ap,Wn*pi*Max_analog_F);
[bz,az] = impinvar(bs,as,Max_analog_F) ;
sys=tf(bz,az,T) ;
[H,W]=freqz(bz,az,512,Max_analog_F);
figure(2)
plot(W,20*log10(abs(H)));
xlabel('f/Hz');ylabel('幅值/dB ');title('滤波器频域幅度谱');
axis([0,ffs,-80,0]);
grid on;
%%
y_noise = 0.3*randn(length(y),2);
noise_y=y+y_noise;
figure(1);
subplot(2,3,2);
plot(Time_scale,noise_y);
title('加噪声时域波形');xlabel('时间长度/秒')
Y_N=fftshift(fft(noise_y));
amp_Y_N=abs(Y_N);