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.
fftfilt()
) B) most of your algorithm makes little sense to me. If would help if you would try to write out the underlying math and assumptions. C) I recommend using proper unit testing and debug techniques: start with simple impulse responses and signals with known answers and work your way. I'm guessing both algorithms are wrong. D) The correct algorithm is indeed overlap-add but that uses no windowing and zero padding which I'm not seeing here. $\endgroup$XFrames
andK
) without zero-padding first to get the correct length. In the time-domain, when you convolve two signals the resulting length is $N_{1} + N_{2} - 1$ where $N_{i}$ is the length (in samples) of each signal vector. You have to “replicate” that in the frequency domain by zero-padding the signal vectors before performing the FFT function. You are not doing that and you end up with circular convolution here, which is incorrect. $\endgroup$