dsp.STFT

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

iscola Function

The iscola function checks that the specified window and overlap satisfy the COLA constraint to ensure that the inverse short-time Fourier transform (ISTFT) results in perfect reconstruction for non-modified spectra. The function returns a logical true if the combination of input parameters is COLA-compliant and a logical false if not. The method argument of the function is set to 'ola' or 'wola' depending on whether the inversion method uses weighted overlap-add (WOLA).

Check if hann() window of length 120 samples and an overlap length of 60 samples is COLA compliant.

winLen = 120;
overlapLen = 60;
win = hann(winLen,'periodic');
tf = iscola(win,overlapLen,'ola')

tf = logical
   1

Initialization

Initialize the dsp.STFT and dsp.ISTFT System objects with this hann window that is COLA compliant. Set the FFT length to equal the window length.

frameLen = winLen-overlapLen;
stf = dsp.STFT('Window',win,'OverlapLength',overlapLen,'FFTLength',winLen);
istf = dsp.ISTFT('Window',win,'OverlapLength',overlapLen,'WeightedOverlapAdd',0);

Reconstruct Data

Compute the STFT of a random signal. Set the length of the input signal to equal the hop length (window length – overlap length). Since the window is COLA compliant, the ISTFT of this non-modified spectra perfectly reconstructs the original time-domain signal.

To confirm, compare the input, x to the reconstructed output, y. Due to the latency introduced by the objects, the reconstructed output is shifted in time compared to the input. Therefore, to compare, take the norm of the difference between the reconstructed output, y and the previous input, xprev. The norm is very small, indicating that the output signal is a perfectly reconstructed version of the input signal.

n = zeros(1,100);
xprev = 0;
for i = 1:100
    x = randn(frameLen,1);
    X = stf(x);
    y = istf(X);
    n(1,i) = norm(y-xprev);
    xprev = x;
end       
max(abs(n))

ans = 
1.4003e-13

ISTFT with Weighted Overlap-Add (WOLA)

In WOLA, a second window called the synthesis window, $\mathit{f}\left(\mathit{n}\right)$, is applied after the IFFT operation and before overlap-add. The synthesis and analysis windows are typically identical and are usually obtained by taking the square root of windows satisfying COLA (thereby ensuring perfect reconstruction).

iscola Function

Check if sqrt(hann()) window of length 120 samples and an overlap length of 60 samples is WOLA compliant. Set the method argument of the iscola function to 'wola'. The output of the iscola function is 1 indicating that this window is WOLA compliant.

winWOLA = sqrt(hann(winLen,'periodic'));
tfWOLA = iscola(winWOLA,overlapLen,'wola')

tfWOLA = logical
   1

Reconstruct Data with WOLA

Release the dsp.STFT and dsp.ISTFT System objects and set the window to sqrt(hann(winLen,'periodic')) window. To use weighted overlap-add on the ISTFT side, set the 'WeightedOverlapAdd' to true.

release(stf);
release(istf);
stf.Window = winWOLA;
istf.Window = winWOLA;
istf.WeightedOverlapAdd = true;

n = zeros(1,100);
xprev = 0;
for i = 1:100
    x = randn(frameLen,1);
    X = stf(x);
    y = istf(X);
    n(1,i) = norm(y-xprev);
    xprev = x;
end       
max(abs(n))

ans = 
4.2068e-15

The norm of the difference between the input signal and the reconstructed signal is very small indicating that the signal has been reconstructed perfectly.