MATLAB带通滤波器：分离歌曲的音频处理

clc
clear
close all
%原始数据导入
[y,fs]=audioread('Song_noisy.wav');%高版本MATLAB的wavread函数已删除,fs为采样频率
t=(0:length(y)-1)/fs;    %时间
figure(1)
subplot(2,1,1)
plot(t,y)
title('原始语音信号的时域波形')
xlabel('秒')
ylabel('幅度')
%原始语音信号FFT
NFFT=length(y);
y1=fft(y,NFFT);
f=(0:NFFT-1)/NFFT*fs;
y1(1,1)=abs(y1(1,1))/2;   %直流量除以2
subplot(2,1,2)
plot(f(1:NFFT/2+1)/1000,abs(y1(1:NFFT/2+1))*2/NFFT)
xlabel('千赫兹')
ylabel('幅度')
title('原始语音信号的幅频谱')
%音频1 范玮琪《最重要的决定》
%巴特沃斯带通滤波器
wc=[2*pi*40000/fs,2*pi*50000/fs];wst=[2*pi*35000/fs,2*pi*55000/fs];
Rp=1;Rst=20;
[N,Wn]=buttord(wc/pi,wst/pi,Rp,Rst);
[B,A]=butter(N,Wn);
[h,w]=freqz(B,A);
figure(2)
subplot(3,1,1)
plot(w/pi,20*log10(abs(h)));
title('巴特沃斯带通滤波器(截取40khz到50khz频率信号)')
xlabel('pi');
ylabel('gain in dB');
y2=filter(B,A,y);%滤波
%相干解调
y2=y2.*cos(2*pi*45e3*t');
%构建低通滤波器
wc=2*pi*2000/fs;wst=2*pi*3000/fs;
Rp=1;Rst=20;
[N,Wn]=buttord(wc/pi,wst/pi,Rp,Rst);
[B,A]=butter(N,Wn);
y2=filter(B,A,y2);%滤波
y2=2*y2;       %理论上相干解调后的信号是1/2*m(t)
subplot(3,1,2)
plot(t,y2)
title('音频1的时域波形')
xlabel('秒')
ylabel('幅度')
%FFT
NFFT=length(y2);
y1=fft(y2,NFFT);
f=(0:NFFT-1)/NFFT*fs;
y1(1,1)=abs(y1(1,1))/2;   %直流量除以2
subplot(3,1,3)
plot(f(1:NFFT/2+1)/1000,abs(y1(1:NFFT/2+1))*2/NFFT)
xlabel('千赫兹')
ylabel('幅度')
title('音频1的幅频谱')
%音频2 陈楚生《天长地久》
%巴特沃斯带通滤波器
wc=[2*pi*60000/fs,2*pi*70000/fs];wst=[2*pi*55000/fs,2*pi*75000/fs];
Rp=1;Rst=20;
[N,Wn]=buttord(wc/pi,wst/pi,Rp,Rst);
[B,A]=butter(N,Wn);
[h,w]=freqz(B,A);
figure(3)
subplot(3,1,1)
plot(w/pi,20*log10(abs(h)));
title('巴特沃斯带阻滤波器(截取60khz到70khz频率信号)')
xlabel('pi');
ylabel('gain in dB');
y3=filter(B,A,y);%滤波.
%相干解调
y3=y3.*cos(2*pi*65e3*t');
%构建低通滤波器
wc=2*pi*1000/fs;wst=2*pi*2000/fs;
Rp=1;Rst=20;
[N,Wn]=buttord(wc/pi,wst/pi,Rp,Rst);
[B,A]=butter(N,Wn);
y3=filter(B,A,y3);%滤波
y3=2*y3;       %理论上相干解调后的信号是1/2*m(t)
subplot(3,1,2)
plot(t,y3)
title('音频2的时域波形')
xlabel('秒')
ylabel('幅度')
%FFT
NFFT=length(y3);
y1=fft(y3,NFFT);
f=(0:NFFT-1)/NFFT*fs;
y1(1,1)=abs(y1(1,1))/2;   %直流量除以2
subplot(3,1,3)
plot(f(1:NFFT/2+1)/1000,abs(y1(1:NFFT/2+1))*2/NFFT)
xlabel('千赫兹')
ylabel('幅度')
title('音频2的幅频谱')
%保存音频
audiowrite('first.wav',y2,fs)%保存音频1
audiowrite('second.wav',y3,fs)%保存音频2