Denoising using stft

[213cy]

新版 matlab 才有有stft 和 istft 两个函数, 这里用 spectrogram 进行 stft, 用 ifft 进行 istft

[213cy]

%% Load sound
[x,Fs] = audioread('spin_ted001.wav');
xn = x(1:37001);

%% STFT parameters
Nfft=1024;
Nwin=768; 
overlap=512; 
win = sqrt(hanning(Nwin,"periodic"));
%iscola(sqrt(hanning(768,"periodic")),512)

%%  Estimate SNR
[Sn,~,~,Pn]= spectrogram(xn,win,overlap,Nfft); 
[Sx,~,~,Px] = spectrogram(x,win,overlap,Nfft); 
% SNR_est=max( Sx.*conj(Sx)./mean( Sn.*conj(Sn) ,2) - 1 ,0); % a posteriori SNR
SNR_est= max( Px./mean(Pn,2) -1 , 0 );
% 1111111111111111111
%alpha=0.5;
%SNR_est=filter((1-alpha),[1 -alpha],SNR_est,[],2); %a priori SNR
% 2222222222222222222
% a = db2mag(-7*[10:-1:0,1:10]);
% SNR_est = conv2(SNR_est,a/sum(a),'same');
% 3333333333333333333
a = db2mag(-7*[10:-1:0,1:10]);
a = a/sum(a);
SNR_est = conv2(a,a,SNR_est,'same');

%%  Compute TF attenuation map 
beta1=0.5;
beta2=0.7;
lambda=1.9; 
mask=max((1-lambda*((1./(SNR_est+1)).^beta1)),0).^beta2;
STFT=mask.*Sx;

%%  Compute inverse STFT and overlapp add
[Ft,Nt] = size(STFT);
Xtemp = ifft([STFT;conj(STFT(Ft-1:-1:2,:))] );
Xs = Xtemp(1:Nwin,:).*win;

hop = Nwin - overlap;
out=zeros((Nt-1)*hop + Nwin,1);
for index=1:Nt 
    ind=(index-1)*hop ;
    out(ind+1:ind+Nwin)= out(ind+1:ind+Nwin)+Xs(:,index);
end

display(['new out ' num2str(randi(999))])

%% method 2
[T2,~,~,P2] = spectrogram(x,win,overlap,Nfft,Fs,'twosided','MinThreshold',-90); 
[T3,~,~,P3] = spectrogram(x,win,overlap,Nfft,Fs,'twosided');
%spectrogram(x,Nwin,overlap,Nfft,Fs,'reassign','MinThreshold',-95);
isequal( T2,T3 ),isequal( P2,P3 ) 

T2(P2==0)=0;
Xtemp2 = ifft(T2);
Xs2 = Xtemp2(1:Nwin,:).*win;

Nt = size(T2,2);
hop = Nwin - overlap;
out2=zeros((Nt-1)*hop + Nwin,1);
for index=1:Nt 
    ind=(index-1)*hop ;
    out2(ind+1:ind+Nwin)= out2(ind+1:ind+Nwin)+Xs2(:,index);
end

return
%% -----------------    Listen results   ------------------------------------
%%
soundsc(xn,Fs);
soundsc(x,Fs);
soundsc(out,Fs);
soundsc(out2,Fs);

%% -----------------    Display Figure   ------------------------------------      
%% show spectrogram
subplot(211)
spectrogram(x,win,overlap,Nfft,Fs); 
%spectrogram(x(9e4:11e4),win,overlap,Nfft,'reassigned');
colorbar(gca,'off')
subplot(212)
spectrogram(out,win,overlap,Nfft,Fs);
colorbar(gca,'off')

[213cy]

https://user-images.githubusercontent.com/10978952/110899136-87d31980-833b-11eb-872e-cdeff6a1579d.mp4

[x,Fs] = audioread('spin_ted001.mp4');

[213cy]

第二个版本:

x=y(min(a):median(a),1);
xn = y([median(a):max(a)],1);
Fs=44100;

%% STFT parameters
Nfft=1024;
Nwin=768; 
overlap=512; 
win = sqrt(hanning(Nwin,"periodic"));

%%  Estimate SNR
[~,~,~,Pn]= spectrogram(xn,win,overlap,Nfft); 
[Sx,~,~,Px] = spectrogram(x,win,overlap,Nfft); 
%%
% A=quantile(Pn,0.35,2).*quantile(Pn,0.65,2)./quantile(Pn,0.2,2)./quantile(Pn,0.8,2);
% noise_A= A.*quantile(Pn,0.5,2);
 noise_B= 2*quantile(Pn,0.35,2);
% noise_C= mean(Pn,2);
%noise_D= median(Pn,2);

noise_r= min( noise_B./Px  , 1 );

fa = db2mag(-3*[20:-1:0,1:20]);
fa = fa/sum(fa);
fb= db2mag(-12*[5:-1:0,1:5]);
fb = fb/sum(fb);
%noise_r = conv2(fb,fa,noise_r,'same');
noise_r = conv2(noise_r,fa,'same');
%noise_r = conv2(noise_r,fb','same');

beta1=0.5;
beta2=0.4;
lambda=logspace(log10(1.1),log10(4.5),Nfft/2+1)'; 
mask=max( 1-lambda.*noise_r.^beta1 , 0).^beta2;
%mask(Nwin/2-30:end,:)=0;
STFT=mask.*Sx;

%%  Compute inverse STFT and overlapp add
[Ft,Nt] = size(STFT);
Xtemp = ifft([STFT;conj(STFT(Ft-1:-1:2,:))] );
Xs = Xtemp(1:Nwin,:).*win;

hop = Nwin - overlap;
out=zeros((Nt-1)*hop + Nwin,1);
for index=1:Nt 
    ind=(index-1)*hop ;
    out(ind+1:ind+Nwin)= out(ind+1:ind+Nwin)+Xs(:,index);
end

display(['new out ' num2str(randi(999))])
return

%% -----------------    Listen results   ------------------------------------
%%
soundsc(xn,Fs);
soundsc(x,Fs);
soundsc(out,Fs);
figure,plot(out)
y=out;
%% -----------------    Display Figure   ------------------------------------      
%% show spectrogram
figure
subplot(211)
spectrogram(x,win,overlap,Nfft,Fs); 
colorbar(gca,'off')
subplot(212)
spectrogram(out,win,overlap,Nfft,Fs);
colorbar(gca,'off')

[213cy]

第5种估计噪音功率的方法

%Pn=rand(5,4);
%%
[r,c]=size(Pn);
Pn_sort=sort(Pn,2);
%plot(Pn_sort')
%
col = ceil(c*0.2);
val =max(Pn_sort(:,col))/(c-col+1)*(c:-1:1);
[~,ind]=min(abs(Pn_sort - val),[],2);
% [~,b]=min(abs(temp - val),[],2,'linear'); matlab2021 Option

F=Pn_sort(ind*r-(r-1:-1:0)');
noise_r= min( F./Px  , 1 );
%%
A=quantile(Pn,0.5,2).*quantile(Pn,0.35,2).*quantile(Pn,0.65,2)./quantile(Pn,0.2,2)./quantile(Pn,0.8,2);
B=2*quantile(Pn,0.5,2);
C=mean(Pn,2);
D=median(Pn,2);
plot([A,B,C,D,F])
legend('Cx','2x','mean','median','linear')