I’m new to DSP programming and I’m trying to learn a bit about convolution and deconvolution and FFTs. I have a project where I am taking a signal (a sine sweep) and convolving it with a shorter Impulse Response. Then I take the convolved signal and deconvolve it with the original signal to recover the IR. I’m writing this in C using the Accelerate Framework’s vDSP functions on OS X. The output from the convolved stage is as expected, though scaled as discussed below, but the output from the deconvolve stage doesn’t match the original IR.
In the convolve stage I get a length for my FFT, N, as close as possible to double the length of the IR like this (because the
vDSP_fft_zrip() function takes a power of 2 parameter for the FFT length):
vDSP_Length log2_N = ceil(log2(irLengthInSamples * 2)); UInt32 N = pow(2.0, log2_N);
I pad the signal out with 0s to be a multiple of N, and pad the IR out with 0s to length N. I iterate through the signal, taking N/2 samples and padding out to length N with 0s, FFT both this padded signal chunk and the padded IR, multiply the results, and IFFT back to the time domain. I scale the results by multiplying by 1/2N to compensate for the scaling introduced by the forward and inverse calls to
vDSP_fft_zrip() . Then I overlap-sum each iteration’s results in to an output array of length signal+IR.
The convolved output signal seems fine, except that it is has a peak loudness of around +10dB (the input signal is at -12dB and the IR has a peak of -15dB). This seems weird, but I don’t know enough to know if this is expected or not. So I normalise the convolved output to 0dB. Normalising or not doesn’t make any significant change to the output of the deconvolve stage..
For the deconvolution stage, I figure the IR length to be the convolved signal length minus the original signal length, and create an FFT length like this:
irLengthInSamples = (convolvedSamplesCount - signalSamplesCount) vDSP_Length log2_N = ceil(log2(rLengthInSamples)); UInt32 N = pow(2.0, log2_N);
I pad the convolved signal out to be a multiple of N, and pad the original signal to the same length. Then I iterate through them in chunks of N length. N samples of the convolved signal and N samples of the original signal are each FFT’d and the frequency values for the convolved signal are divided by the frequency values of the original signal. The results are IFFT’d, scaled again by 1/2N and summed into an output array of length N.
But the output does not match the original IR. The image below shows the original IR waveform on top and the deconvolution output on the bottom.
One issue I'm aware of is that because the original signal is padded at the end with 0s, the last iteration through the deconvolve loop includes division by 0, producing NaNs. I’ve been handling this by replacing NaNs with the equivalent samples from the convolved signal before summing them, but I suspect I shouldn’t be getting into this situation in the first place.
Can anyone point out where I'm going wrong? As I said, I'm new to DSP programming so I suspect I may be making a fundamental mistake somewhere.