[x,fs11]=audioread('I:\舊數據\xinhao\ke she\第五\shi miao.mp3');
fs=1000;
%sound(x,fs11); %聲音
%x=x(:,1); %單聲道
Fs=fs; %頻率
N=length(x); %音樂的長度
n=0:N-1;
f=n*fs/N; %頻率序列
tt=(N-1)/fs;
ts=1/fs;
t=0:ts:tt;
X=fft(x,N);
figure(1);
subplot(211);plot(t,x);title('音頻時域波形');
xlabel('t');ylabel('x');
subplot(212);plot(f,X);title('音頻頻域波形');
xlabel('f/hz');ylabel('X');
%設計噪聲
nt1=cos(2*pi*300*n/fs); %第一種噪聲
%nt1=nt.';
nt2=2*randn(1,N)-1; %第二種噪聲
nt=nt1+nt2;
nt=nt.';
%fp=23000;f_s=18000;
fp=400;f_s=390;
wp=fp*2/fs; ws=f_s*2/fs
Ap=3;As=18;
[N1,wn]=buttord(wp,ws,Ap,As);
[b,a]=butter(N1,wn,'high');
[h,w]=freqz(b,a);
ntt=filter(b,a,nt);
%ntt=ntt.';
Ntt=fft(ntt,N);
figure(2);
plot(f,Ntt);title('噪聲頻域波形');
xlabel('f/hz');ylabel('Ntt');
%加噪過程
x_ntt=x + ntt; %音頻+噪聲
%x_ntt=x_ntt.';
%sound(x_ntt,fs11);
x_ntt_fft=fft(x_ntt,N);
figure(3);
subplot(211);plot(t,x_ntt);title('噪聲時域波形');
xlabel('t');ylabel('x-noise');
subplot(212);plot(f,x_ntt_fft);title('加噪音頻頻域波形');
xlabel('f/hz');ylabel('x-noise-fft');
%設置濾波器
Ap=3;As=18;
fp1=300;f_s1=350;
%fp1=15000;f_s1=17000; %觀察加噪音頻頻帶和音頻頻帶得到濾波器阻帶通帶的截止頻率
wp1=fp1*2/fs; ws1=f_s1*2/fs
Ap=3;As=18;
%[Nn,w2]=cheb1ord(wp1,ws1,Ap,As); %切比雪夫濾波器
%[b1,a1]=cheby1(Nn,Ap,w2);
%[Nn,w2]=ellipord(wp1,ws1,Ap,As); %橢圓濾波器
%[b1,a1]=ellip(Nn,Ap,As,w2);
[Nn,w2]=buttord(wp1,ws1,Ap,As); %巴特沃斯濾波器
[b1,a1]=butter(Nn,w2);
[H1,w3]=freqz(b1,a1);
figure(4);
plot(w/pi,abs(H1));title('濾波器幅頻圖');
xlabel('w/pi');ylabel('H1');
%進行濾波
x2=filter(b1,a1,x_ntt);
N2=length(x2);
n2=0:N2;
f2=n*fs/N2;
X2=fft(x2,N2);
sound(x2,fs11);
figure(5);
subplot(211);plot(t,x2);title('濾波后的音頻時域圖');
xlabel('f/hz');ylabel('x');
subplot(212);plot(f2,X2);title('濾波后的音頻頻域圖');
xlabel('f/hz');ylabel('X');
|