Why will the output be imprinted with delay effect or echo after convolving two direct dry signals' FFT values (convolution filtering)

Question

I am recently working on the reproduction of the filtering effect of target play back devices, like phone, or speakers, etc using convolution techniques in MATLAB.

I firstly created a function called "convolutionFilter.m", which essentially performs a standard convolution operation by multiplying two FFT values

function [z,x,Fs]= convolutionFilter(inputfileName,kernelfileName, N, overlap)

% excerpt the first 4096 samples long's segment as the filter kernel, whcih
% is through 1~ N.  
  [k, Fs] = audioread(kernelfileName, [1, N]); 

% make the kernel a mono signal
  k = k(:,1);

% apply hann window before all the FFT algorithm, in order to get more 
% realistic results and reduce spectral leakage
  hw = hann(N);   

% get the FFT results of the filter kernel
   K = fft(k .* hw);  

% load an input signal
   [x, Fs] = audioread(inputfileName);
 
% make the input sigal a mono signal as well
   x = x(:,1);

% break the input signal into subframes, make sure N is as long as the
% kernel "k" 
   xFrames = audioFrames(x, N, overlap);

% take the FFT result of all windowed subframes, or say,get the FFT result 
% of xFrames by adding each subframe's FFT 
   XFrames = fft(xFrames .* hw);

% multiply the FFT of xFrames by the FFT of filter kernel to get the
% frequency domain, which is the step of convolution
   ZFrames = K .* XFrames; 
   
% get the output signal by using IFFT and real function from its frequency 
% spectrum back to the time domain, which will be the audio samples with 
% all the real results  
   zFrames = real(ifft(ZFrames));  
   
% since we have got this new filtered value, we need to apply hann window 
% to it
   zFrames = zFrames .* hw; 

% use overlap-add method to construct the filtered output signal 
   z = frameAssembler(zFrames, overlap); 

% normalize the output signal
   z = z/max(abs(z));  

% make sure the new filtered output signal and the input signal 
% are the same magnitude
  z = z * max(abs(x));  
  
end

And the output signal returned by the above function will bring some delay or echo, which is unwanted. So, I created another function which is not the multiplication between the two's FFT values, but firstly get the magnitude response of the Impulse response(here is called kernel), and then get its filtered array index between 0 to 1 by using the kernelMagnitude array to divide the maximum single value amongst the kernelMagnitude array, followed by the final multiplication between the FFT of input signal and this filterCurve or say energy scaler.

This is the second function

function [z,x,Fs] = convolutionFiltercurve(inputfileName,kernelfileName, N, overlap)
   
% cut the first 4096 samples long's segment as the filter kernel,from 1~ N  
  [k, Fs] = audioread(kernelfileName, [1, N]); 
  
% make the kernel a mono signal
  k = k(:,1);

% apply hann window before all the FFT algorithm, in order to get more 
% realistic results and reduce spectral leakage
   hw = hann(N);   

% get the FFT result of the filter kernel
   K = fft(k .* hw);  

% load an input signal
   [x, Fs] = audioread(inputfileName);
 
% make the input signal a mono signal
   x = x(:,1);

 % break the input signal into subframes, make sure N is as long as the
 % kernel "k" 
  xFrames = audioFrames(x, N, overlap);

% take the FFT result of all windowed subframes, or say,get the FFT result 
% of xFrames by adding each subframe's FFT value 
   XFrames = fft(xFrames .* hw);

% get the kernel's magnitude spectrum which is also the absolute FFT 
% result's array of "K"
   kernelMagnitude = abs(K);

% normalize the kernelMagnitude array to get an array of numbers between 
% 0 to 1 by using the kernelMagnitude array to divide the maximum single
% value of the kernelManitudet array.Filnally, we can get an array of 
% proportion through 0 to 1, which consists of the filterCurve we want
 filterCurve = kernelMagnitude/max(kernelMagnitude); 

% multiply this filterCurve array by the XFrames FFT array to get the final
% filterd Frames' frequency spectrum(FFT),because at this moment, the 
% filterCure is a magnitudeScalar
 ZFrames = filterCurve .* XFrames;
   
% get the output signal by using IFFT and real function from its frequency 
% spectrum back to the time domain, which will be the audio samples with 
% all the real results  
   zFrames = real(ifft(ZFrames));   
   
% since we have got this new value, we need to apply hann window to it
  zFrames = zFrames .* hw;  

% use overlap-add method to construct the filtered output signal  
  z = frameAssembler(zFrames, overlap);  

% normalize the output signal
  z = z/max(abs(z));  

% make sure the new filtered output signal and the input signal 
% are the same magnitude
  z = z * max(abs(x));   
  
end

And the output returned by the second method" convolutionFilterCurve.m" has no delay or echo.

Note: Both methods use the same window size and overlap size when executing them, and both the recorded impulse response (kernel) and the input signal are free from delay or echo. Therefore, I'm curious about the specific reason for the output returned by the first method causing delay, while the output returned by the second method has no delay effect at all? Thanks in advance!

Lastly, the "xFrames = audioFrames(x, N, overlap)" and "z = frameAssembler(zFrames, overlap)" are another two auxiliary functions which correctly perform the overlap-add decomposition and overlap-add reconstruction.

Answer 1

The difference in the behavior of the two methods, particularly the delay or echo effect in the first method, is primarily due to the way the kernel is being used in both methods.

In your first method, you are using the FFT of the entire kernel to convolve with your input signal. If the impulse response (or kernel) has a large number of samples (4096 samples in your case), it means that there is a large time span that is being used to convolve with your input signal. In simple terms, the convolution operation in the time domain is equivalent to applying a delay and weighted sum of the input signal samples. If the kernel has a large number of non-zero samples, it will create a delayed version of the input signal which can sometimes make it sound like an echo.

In your second method, you are using the magnitude spectrum of the kernel (or the absolute FFT result), and then normalizing it to be between 0 and 1. This step essentially flattens the phase spectrum of the kernel and only keeps the magnitude information. The phase of a signal is what gives the signal its time-dependent characteristics such as delay. By removing or flattening the phase information, you are getting rid of the delay characteristics of the kernel, hence the output of the second method doesn't have any delay or echoing effect.

It's important to remember that while the second method might not have the delay or echo effect, it is not accurately reproducing the filtering effect of the target playback devices since the phase information, which is crucial for accurate reproduction of the audio characteristics of the device, is being discarded. Although, it is true that the human ear is not sensitive to phase, this does not mean that we should discard it entirely.

Answer 2

The correct way of doing this is overlap add. Neither version of your code looks remotely like a correct overlap add implementation, so chances both are wrong.

The proper way to determine this is to compare against a "known-good" reference implementation and/or start with simple signals and impulse response where the answer is known and/or an be easily calculated by hand. That's also useful in debugging your code and figuring out where exactly the problems are. Below is code example that does that.

%% script to demonstrate overlap add
nx = 48000; % number of samples
nh = 1024; % impulse response length
x = randn(nx,1);
h = randn(nh,1);
reference = fftfilt(h,x); % reference
y = myOLA(x,h);
d = y-reference;
fprintf("error is %6.2fdB\n",10*log10(sum(d.^2)./sum(reference.^2)));

function y = myOLA(h,x)
nx = length(x);
nh = length(h);
nFFT = nh; % in practice you would choose a "good" size here
numFrames = ceil(nx/nFFT);  % number of frames
y = zeros(numFrames*nFFT,1);  % init output
overlap = zeros(nFFT,1);      % init overlap
% main loop over all frames
cnt = 0;
fh = fft(h,2*nFFT); % filter spectrum zero pad and FFT
t = 1:nFFT; % frame time vector
for i = 1:numFrames
  if i == numFrames  % handle the last frame which may be short
    x1 = x(cnt+1:end,:);
  else
    x1 = x(cnt+t,:); % get frame
  end
  fx1 = fft(x1,2*nFFT);  % zero pad and FFT
  fy1 = fx1.*fh;  % mulitply spectra
  y1 = real(ifft(fy1)); % invers FFT
  y(cnt+t,:) = y1(t,:) + overlap; % create output
  overlap = y1(nFFT+t,:); % save overlap for next frame
  cnt = cnt + nFFT;  % update position index
end % for
end % function