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 (
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).
You can download the input audio files at these links
I wrote two OCTAVE procedures
- with direct convolution using
- 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).
I noticed that the result of STFT procedure is quite different from direnct convolution/FFT product procedure. You can download the output file at this link http://www.sun-art.org/download/matlab/STFTAudioConv.wav
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'); [x2, FS] = audioread('conv_test\test_02.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, FS] = audioread('conv_test\test_01.aif'); [x2, FS] = audioread('conv_test\test_02.aif'); % 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
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
- take 1024 (L) samples of inputs with windowing
- zero padding with (N-L) = 1024 zeros. The buffer length in the FFT is now 2048 samples long.
- FFT of zero padded inputs
- FFT product
- overlapp add the output
- HopSize advance
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, FS] = audioread('conv_test\test_01.aif'); [x2, FS] = audioread('conv_test\test_02.aif'); % 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 % ================== OVERLAPP ADD ================== 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