%% Short-Time Spectral Attenuation
% Short-time spectral attenuation is achieved by applying a time-varying attenuation to the short-time spectrum of a noisy signal. The gain of the attenuation is determined by the estimate of the noise power in each subband of the spectrum. This gain, when applied to the noisy spectrum, attenuates the subbands with higher noise power and lifters the subbands with lesser noise power.
% Here are the steps involved in performing the short-time spectral attenuation:
% Analyze the noisy input signal by computing the short-time Fourier transform (STFT).
% Multiply each subband of the transformed signal with a real positive gain less than 1.
% Synthesize the denoised subbands by taking the inverse short-time Fourier transform (ISTFT). The resconstructed signal is the denoised input signal.
% Use the  dsp.STFT and dsp.ISTFT objects to compute the short-time and the inverse short-time Fourier transforms, respectively.
%
%% Noisy Input Signal
% The input is an audio signal sampled at the 22,050 Hz. The dsp.AudioFileReader object reads this signal in frames of 512 samples. The audio signal is corrupted by white Gaussian noise that has a standard deviation of 0.05. Use the audioDeviceWriter object to play the noisy audio signal to your computer's audio device.

FrameLength = 512;
afr = dsp.AudioFileReader('speech_dft.wav',...
    'SamplesPerFrame',FrameLength);
adw = audioDeviceWriter('SampleRate',afr.SampleRate);

noiseStd = 0.05;
while ~isDone(afr)
    cleanAudio = afr();
    noisyAudio = cleanAudio + noiseStd * randn(FrameLength,1);
    adw(noisyAudio);
end
reset(afr)

%
%% Initialize Short-Time and Inverse Short-Time Fourier Transform Objects
% Initialize the dsp.STFT and dsp.ISTFT objects. Set the window length equal to the input frame length and the hop length to 16. The overlap length is the difference between the window length and the hop length, OL = WL – HL. Set the FFT length to 1024.

WindowLength = FrameLength;
HopLength = 16;
numHopsPerFrame = FrameLength / 16;
FFTLength = 1024;

% The window used to compute the STFT and ISTFT is a periodic hamming window with length 512. The ConjugateSymmetricInput flag of the istf object is set to true, indicating that the output of the istf object is a conjugate-symmetric signal.

win = hamming(WindowLength,'periodic');
stf = dsp.STFT(win,WindowLength-HopLength,FFTLength);
istf = dsp.ISTFT(win,WindowLength-HopLength,1,0);

%
%% Gain Estimator
% The next step is to define the gain estimator parameters. This gain is applied to the noisy spectrum to attenuate the subbands with higher noise power and lifter the subbands with lesser noise power.

dec = 16;
alpha = 15;
stftNorm = (sum(win.*win) / dec).^2;

%
%% Spectral Attenuation
% Feed the audio signal to stf one hop-length at a time. Apply the estimated gain to the transformed signal. Reconstruct the denoised version of the original speech signal by performing an inverse Fourier transform on the individual frequency bands. Play the denoised audio signal to the computer's audio device.

while ~isDone(afr)
    cleanAudio =  afr();
    noisyAudio = cleanAudio + noiseStd * randn(FrameLength,1);
    y = zeros(FrameLength,1); % y holds the denoised audio frame
    
    % Feed audio to stft one hop-length at a time
    for index = 1:numHopsPerFrame        
        X = stf(noisyAudio((index-1)*HopLength+1:index*HopLength));        
        % Gain estimator
        Z = abs(X).^2 / (noiseStd^2 * alpha) / stftNorm;
        Z(Z<=1) = 1;
        Z = 1 - 1./Z;
        Z = sign(Z) .* sqrt(abs(Z));
        X = X .* Z;        
        % Convert back to time-domain
        y((index-1)*HopLength+1:index*HopLength) = istf(X);        
    end    
    % Listen to denoised audio:
    adw(y);
end

% Copyright 2018 The MathWorks, Inc.