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:
- 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