[PD] deconvolution in pd

[Charles Henry]

I have implemented the fourier domain adaptation algorithm and worked
out most of the bugs.  It is doing a good job of restoring frequency
and phase response for short (8192) length filtering on blocks of
16384 samples.  The two signals, one reference signal and the received
signal from the microphone line up on top of one another pretty well.
The problem, which you can see in the screenshot, is that the signal
at the beginning and end are getting left out.  This may have
something to do with my handling of the frequencies at zero and at the
Nyquist frequency (of the wavelet subband being adapted).  The FT at
these two frequencies was being problematic, so instead of dividing by
the numbers (which sometimes turn out as 0) I used different lines:
 mag[0]=((fabsf(refsub[0]) > fabsf(recsub[0]))? 1.1 : 0.9);
 mag[half]=((fabsf(refsub[half]) > fabsf(recsub[half]))? 1.1 : 0.9);
these circumvent the problem with division by zero, and tend to
converge...but as you can see, they don't coverge to the right number.
 Any ideas what I could replace these lines with?

I'll enclose the code for the adaptive filters called controlfreak~.c
also, and see if any suggestions for improving the code come up (this
is a very high latency filtering routine; once I get it fully working
with high latency, I'll try to modify it to reduce the latency)

A representative screenshot can be found here:
http://czh.zapto.org/public/adaptation.jpg

Algorithm:
> my idea was this: given alpha <=1 (an adaptation coefficient)
> we do:
> 1. Fourier transform of 3 things: reference signal (length 2N),
> received signal from microphone (length 2N), and the weights (length
> N, padded with N zeros at the end).
> 2. divide (FT of reference) by (FT of received)
> 3. unwrap phase; take result of (2) as phase and magnitude
> 4. multiply (FT of weights) by mag^alpha*( cos (phase*alpha) + i
> sin(phase*alpha) )
> 5. take IFT of weights, and throw away final N samples
Lastly, the changes to this:
6. multiply the weights by a gaussian centered at the average time
location, E(t), handled in a periodic sense; the standard deviation is
a function of alpha, so that the gaussian converges to 1, everywhere,
as alpha goes to 0
7. measure and renormalize weights
8. measure and reject adaptation if correlation decreases; if
correlation improves, but is negative, multiply weights by -1

Chuck