语音降噪算法：基于短时傅里叶变换的谱减法

addpath('./STFT')
clc
clear all
close all

[x,~]=audioread('C:/Users/ASUS/Desktop/20230607_210020.wav');          %读取纯净语音信号
[d,fs]=audioread('C:/Users/ASUS/Desktop/20230607_210020.wav');       %读取噪声信号
x=x(:,1);
d=d(:,1);
nfft=256;
stepsize=0.5;
snr=5;

d=sqrt(norm(x)^2/(10^(snr/10)*norm(d)^2))*d; %考虑信噪比
y=[zeros(fs*1,1);x(1:fs*9)]+d(1:fs*10);  %生成带噪信号，前1s 为纯噪声段

Y=stft(x,nfft,stepsize*nfft,1);%带噪信号的时频谱
Ya=abs(Y);%幅度谱
Yp=angle(Y);%相位谱
Xa=zeros(size(Ya));

%带噪语音谱平滑
Ys=Ya;
for i=2:(size(Ya,2)-1)
    Ys(:,i)=(Ya(:,i-1)+Ya(:,i)+Ya(:,i+1))/3;
end

N=mean(Ya(:,1:10),2); %噪声谱初始化 使用前0.32s的纯噪声段
N_max=zeros(size(N));  %最大噪声残留初始化
alpha=0.9;%噪声谱更新 平滑系数
beta=0.1;

%静音段噪声谱估计和更新
for i=1:floor(fs/256)
    N=alpha*N+(1-alpha)*Ya(:,i);
    N_max=max(N_max,Ya(:,i)-N);
    Xa(:,i)=beta*Ya(:,i);
end

%带噪语音段降噪
for i=(floor(fs/256)+1):size(Y,2)-1
     X1=Ys(:,i)-N;  %语音谱的初始估计    
	 for j =1:length(X1)
	     if X1(j)<N_max(j)
	         X1(j)=min(X1(j),min(Ys(j,i-1)-N(j),Ys(j,i+1)-N(j)));
	     end
	 end
    Xa(:,i)=max(X1,0);
end

Xa(:,end)=max(Ys(:,end)-N,0);

%反变换得到增强语音
X=Xa.*exp(1i*Yp);
x_est=istft(X,nfft,stepsize*nfft,1);
x_est=x_est(:);

figure
subplot(311)
plot((0:length(x)-1)/fs,[zeros(fs,1);x(1:fs*9)]);
title('纯净语音')
subplot(312)
plot((0:length(y)-1)/fs,y);
title('带噪语音')
subplot(313)
plot((0:length(x_est)-1)/fs,x_est);
title('增强语音')
figure
subplot(311)
spectrogram(x(1:fs*9),hann(256),50,256,fs,'yaxis')
title('纯净语谱图')
subplot(312)
spectrogram(y(fs+1:end),hann(256),50,256,fs,'yaxis')
title('带噪语谱图')
subplot(313)
spectrogram(x_est(fs+1:end),hann(256),50,256,fs,'yaxis')
title('增强语谱图')