# What's wrong with my STFT implementation?

I'm studying the linear convolution for achieve more complex results in Audio FX.

My main goal is to implement a realtime convolution between two stream of sound in live perfomance (for example Sound 1 from a synthesizer and Sound2 from a guitar o microphone) .

The length of sounds are not determinate and for testing i recorded two part of these streams. You can think the file used in Matlab as a part of infinite streams.

I have two audio files (x1 and x2) with same fixed length (len1 and len2 = 673792 samples) and from DSP theory i know that the length of convolution array is lenOut = (len1+len2-1).

I wrote two OCTAVE procedures

1. with direct convolution using conv(n,m) function
2. with FFT product of input signals with zero padding

Both procedures gives the same results. You can download the output audio files (30 seconds long each) at these links:

The next step is to make a realtime convolution using STFT but before going in C++ i would like to improve it and test it in Matlab/Octave. (my reference are: "SPECTRAL AUDIO SIGNAL PROCESSING" by Julius O. Smith III and "DAFX: Digital Audio Effects" by Udo Zölzer).

The file is 30 seconds long, but from 15 seconds onwards the audio is zero. It'possible to achieve the result of function Matlab conv() with STFT procedure ? What's wrong with my STFT ?

Below you will find the three procedures.

1) LINEAR CONVOLUTION

    [x1, FS] = audioread('conv_test\test_01.aif');

y        = conv(x1,x2);
outFile  = y / max(abs(y));

audiowrite ('conv_test\conv.wav', outFile, FS );  %Octave


2) FFT PRODUCT

    % ============= LOAD audio files x1 and x2==============

% x1 and x2 are files with the same length / length(x1) = length(x2) = 673792 samples
len1   = length(x1);
len2   = length(x2);
lenOut = len1 + len2 - 1;

N      = 2^nextpow2(lenOut);   %for FFT and padding

% ============= PADDING _ & scaling =============
inx1   = [x1;  zeros(N-len1,1)] / max( abs(x1) );
inx2   = [x2;  zeros(N-len2,1)] / max( abs(x2) );

%============= FFT of signal x1 =============

f1     = fft(inx1);
r1     = abs(f1);
phi1   = angle(f1);
fft1   = (r1.* exp(i*phi1));

%============= FFT of signal x2 =============

f2     = fft(inx2);
r2     = abs(f2);
phi2   = angle(f2);
fft2   = (r2.* exp(i*phi2));

%============= FFT PRODUCT  and iFFT =============

Yfft   = fft1 .* fft2;
yout   = real( ifft(Yfft) );

%============= Write and Save iFFT (linear convolution) =============

outFile = yout(1:lenOut) / max( abs(yout) );
audiowrite ('conv_test\convFFT.wav', outFile, FS);  %Octave


3) STFT

For Overlapp-Add i define an Hanning window with L = 1024 samples, an Hopsize of L/2 = 512 samples and an FFT size N = 2^nextpow2(L+1) = 2048 samples. In the while loop

1. take 1024 (L) samples of inputs with windowing
2. zero padding with (N-L) = 1024 zeros. The buffer length in the FFT is now 2048 samples long.
3. FFT of zero padded inputs
4. FFT product
5. iFFT

Fixing the windows lenght L = 1024 implicitly I also fixed the signal (x1) and filter (x2) size L and FFT size = 2*L .

My inpulse response (x2) as the other signal is varying-time. My processor has two input buffer (left and right channel) for example 16 samples long. I collect L = 1024 sample in new two auxilaty buffers (in1 and in2 in Script 3) and after windowing i extend them with zero-pad until the FFT size ( = 2*L)

    % ============= LOAD audio files x1 and x2==============

% x1 and x2 are files with the same length / length(x1) = length(x2) = 673792 samples
len1     = length(x1);
len2     = length(x2);
lenOut   = len1 + len2 - 1;

L         = 1024;               % length of analisys windows
w1        = hanningz(L);        % w = .5*(1 - cos(2*pi*(0:L-1)'/(L)));
hop       = L/2;
N         = 2^nextpow2(L+1);    % FFT length

out       = zeros(lenOut,1);    %init of output
start     = 0;                  %index of while loop
pout      = 0;
pend      = (length(x1) - L);   %end of while loop

while start<pend

in1    = x1(start+1:start+L) .* w1;       %windowing
in2    = x2(start+1:start+L) .* w1;

max1   = max(abs(in1));
max2   = max(abs(in2));

buf1   = [in1; zeros(N-L,1)] / max1;  % zero padding
buf2   = [in2; zeros(N-L,1)] / max2;

f1     = fft(buf1);                   % FFT input 1
r1     = abs(f1);
phi1   = angle(f1);
ft1    = (r1.* exp(i*phi1));

f2     = fft(buf2);                   % FFT input 2
r2     = abs(f2);
phi2   = angle(f2);
ft2    = (r2.* exp(i*phi2));

YConv  = ft1 .* ft2;                  % Convolution - FFT Product

y      = real( ifft(YConv) );         % iFFT

out(pout+1:pout+N) = out(pout+1:pout+N) + y;  % Overlapp-Ddd

start  = start + hop;
pout   = pout + hop;

end

% ================== Write and Save ==================

out = out / max(abs(out));
audiowrite ('conv_test\STFTAudioConv.wav', out, FS ); %Octave

• You are mixing up STFT and Overlap-Add. STFT is used to visualize Amplitude v/s Frequency v/s Time. It is different from obtaining convolution of 2 signals using Overlap-Add. You are splitting second sequence also. That is not correct. Overlap - Add is used when one sequence is large length and another small length. You do not split second sequence. But here both seems to be having large lengths. Apr 9 '20 at 10:26
• @jithin you can certainly also split the second sequence, if the impulse response is very long. That's a very common thing to do with things like room impulse responses. It's typically called a block convolver. Apr 9 '20 at 13:35
• Thanks @Hilmar. Never heard of it before. Will explore. Apr 9 '20 at 13:41
• I would suggest that the question is re-phrased to something a bit more specific. Perhaps it could become about clarifying the methods at a higher level? This will make you see if something you are doing in your code should be modified.
– A_A
Apr 9 '20 at 15:19
• @A_A read the second Answers :) Apr 9 '20 at 15:39

A few points

1. It is very unusual to convolve two audio files. You typically convolve one audio file with the impulse response of a filter. Thinking about "long file" and "short-ish" impulse response makes it easier to understand
2. For overlap add the FFT size is determined by the length of the filter. For example, if your file length is 670k samples and your impulse response is 16k samples, you would chose an FFT size of twice the short one, i.e. 32k. You multiply the FFT of each input frame with the FFT of the entire filter.
3. You don't window. Or to be precise: you use a rectangular window.
4. If your impulse response is too long for a single FFT (for latency, memory or CPU reasonse), you can chop this up as well, but that's more complicated and you should not do this until you have Overall Add figured out.
5. Overlap add is NOT the same as STFT. (see jithin's comment)
6. Your code has a lot of unnecessary stuff in it. You extract amplitude and phase from the FFT of the block and put them together again. Why? $$ft1 = f1$$
7. You normalize the input to the maxium of a windowed frame. That's just plain wrong