常见噪声仿真

clear;
clc;
%产生带通噪声（用窄带高斯白噪声做模拟）
%初始条件设置
fs=44100;           %采样频率：44.1kHz
f_howling=10000;    %单频啸叫的对应频率
duration=1;         %音频时长1s
%保存的语音信号路径及名称，助教学姐验证的时候注意一下路径
filename='D:\\project_record\\speech.wav';
filename_lbn='D:\\project_noise\\speech_lbn.wav';
filename_hbn='D:\\project_noise\\speech_hbn.wav';

%读取原始音频信号
speech_record=audioread(filename);
video_len=length(speech_record);

%产生高斯白噪声
noise_wgnl=wgn(fs*duration,1,10,0.06);
noise_wgnh=wgn(fs*duration,1,1,0.11);

%将高斯白噪声通过低通滤波器得到低频噪声
Wp_l1=(20*2*pi)/fs;Wp_h1=(360*2*pi)/fs;
Ws_l1=(30*2*pi)/fs;Ws_h1=(350*2*pi)/fs;

Wc1=[(Wp_l1+Ws_l1)/2/pi,(Wp_h1+Ws_h1)/2/pi];% 截止频率在通带和阻带边界频率的中点
%Wc=(Wp+Ws)/2;
% 计算过渡带宽度
Wb1=Ws_l1-Wp_l1;

%blackman窗设计
M=ceil(5.56*pi/Wb1);
N=2*M;
numl=fir1(N,Wc1,'band',blackman(N+1));
noise_lbn=filter(numl,1,noise_wgnl);
speech_lbn=speech_record+noise_lbn;
%求取信噪比
snr_lbn=snr(speech_record,noise_lbn);
disp('加入低频带通噪声信号的信噪比(dB)：'),disp(snr_lbn)

%将高斯白噪声通过高通滤波器得到高频噪声
Wp_l2=(10000*2*pi)/fs;Wp_h2=(11050*2*pi)/fs;
Ws_l2=(10050*2*pi)/fs;Ws_h2=(11000*2*pi)/fs;

Wc2=[(Wp_l2+Ws_l2)/2/pi,(Wp_h2+Ws_h2)/2/pi];% 截止频率在通带和阻带边界频率的中点
%Wc=(Wp+Ws)/2;
% 计算过渡带宽度
Wb2=Ws_l2-Wp_l2;

%blackman窗设计
M=ceil(5.56*pi/Wb2);
N=2*M;
numh=fir1(N,Wc2,'band',blackman(N+1));
noise_hbn=filter(numh,1,noise_wgnh);
speech_hbn=speech_record+noise_hbn;
snr_hbn=snr(speech_record,noise_hbn);
disp('加入高频带通噪声信号的信噪比(dB)：'),disp(snr_hbn)

%保存录音
audiowrite(filename_lbn,speech_lbn,fs);
audiowrite(filename_hbn,speech_hbn,fs);
%听录音
%sound(speech_lbn,fs);
%sound(speech_hbn,fs);


%绘制图像
%低频带通噪声时域与功率谱图像
t=(0:1:video_len-1)/video_len*duration;%时间轴单位s
n=0:video_len-1;        
N=video_len;      
figure(1)
subplot(2,1,1);
plot(t,noise_lbn);
title('低频带通噪声时域图像');
xlabel('时间：s');
grid on;
R=xcorr(noise_lbn,noise_lbn);
Y=fft(R);
subplot(2,1,2);
plot(10*log10(abs(Y)));
title('低频带通噪声音频信号功率谱图像');

%粉红噪声时域与功率谱图像
video_len=length(speech_record);
t=(0:1:video_len-1)/video_len*duration;%时间轴单位s
n=0:video_len-1;        
N=video_len;      
figure(2)
subplot(2,1,1);
plot(t,noise_hbn);
title('高频带通噪声时域图像');
xlabel('时间：s');
grid on;
R=xcorr(noise_hbn,noise_hbn);
Y=fft(R);
subplot(2,1,2);
plot(10*log10(abs(Y)));
title('高频带通噪声音频信号功率谱图像');



%绘制原始音频信号的时域频域图像
t=(0:1:video_len-1)/video_len*duration;%时间轴单位s
n=0:video_len-1;        
N=video_len;      
%=======时域图像======
figure(3)
subplot(2,1,1);
plot(t,speech_record);
title('原始音频信号时域图像');
xlabel('时间：s');
grid on;
%=======频域图像======     
Y1=fft(speech_record,N);    
mag=abs(Y1);            
f=n*fs/N;               
%取1/2作图
subplot(2,1,2);
plot(f(1:fix(N/2)),mag(1:fix(N/2)));
title('原始音频信号频谱图');
xlabel('频率/Hz');      
ylabel('幅度');         
grid on;    


figure(4)
%绘制加入低频带通噪声的音频信号的时域功率谱图像
%=======时域图像======
subplot(2,1,1);
plot(t,speech_lbn);
title('加入低频带通噪声音频信号时域图像');
xlabel('时间：s');
grid on;
%=======功率谱图像======     
%取1/2作图
subplot(2,1,2);
R=xcorr(speech_lbn,speech_lbn);
Y=fft(R);
subplot(2,1,2);
plot(10*log10(abs(Y)));
title('加入低频带通噪声音频信号功率谱图像');


%绘制加入粉红噪声的音频信号的时域功率谱图像
figure(5)
%=======时域图像======
subplot(2,1,1);
plot(t,speech_hbn);
title('加入高频带通噪声音频信号时域图像');
xlabel('时间/s');
grid on;
%=======功率谱图像======                
R=xcorr(speech_hbn,speech_hbn);
Y=fft(R);
subplot(2,1,2);
plot(10*log10(abs(Y)));
title('加入高频带通噪声音频信号功率谱图像');