0
$\begingroup$

I am developing C++ code to do some modeling that will take as input a parametric model of the complex frequency response (amplitude, phase) and produce as output an impulse response, e.g. using the iFFT. I would like to compare this to a "windowed" impulse response that I have measured for a loudspeaker in a room. The measured impulse can be characterized by an initial portion that corresponds to sound traveling directly through the air to the microphone that lasts a few millisecond, after which sounds that reflect off of the surfaces in the room (floor, walls, etc.) contaminates the signal. The uncontaminated portion is only 3% of the total measured impulse. For instance if I record 16k points at 96kHz, its only the first 470 points of the impulse. I want to compare the impulse response that the iFFT generates from the frequency response model to ONLY the uncontaminated portion of the measured impulse.

I could calculate the entire impulse response from the frequency response and then just throw away 97% of the result but this seems to be very inefficient. The iFFT will be calculated many, many times (thousands probably) while my model is being optimized so I want to make sure that I can make it as efficient as possible. At this point my only option seems to be using FFTW and then just throwing away the data that is not needed (for lack of a better idea).

Is there a fast way to calculate the inverse FFT only for those 470 time points of interest, e.g. not for the entire time span that the FFT can access? I am not intimately familiar with the computation of the FFT and iFFT, so I don't have insight on the answer to this question.

Improve this question
$\endgroup$
0
$\begingroup$

If your inverse FFT is producing a time series that is N times longer than you need, downsample your frequency domain signal by N and then take the IFFT. Try it!

Edit: As has been pointed out, this method will produce undesired time aliasing unless the content after the window of interest is very small or non-existent. This method may still be favorable in cases where the system response decays exponentially. However, if your window ends before the response has decayed by say, 20 or 30 dB, you should consider another method that will be free of aliasing. Either that, or scale your window down only to the point where you are not negatively affected by aliasing, and then throw out any un-needed samples from there.

Improve this answer
$\endgroup$
  • $\begingroup$ Now why didn't I think of that???! Great suggestion. Do I need to be thinking about aliasing (or something like that) when doing the inverse transform? For instance, do I need to make sure that I keep the number of samples high enough so that the impulse dies out within the time domain produced by the iFFT? Are there any other considerations that I should keep in mind when choosing the number of frequency domain samples? $\endgroup$ – Charlie Oct 3 '13 at 23:22
  • 2
    $\begingroup$ @Charlie: Yes, you will have time-domain aliasing if you downsample before the inverse DFT. If the other unwanted samples have negligible values, then you can get away with this approach. However, if they have significant amplitude, then by downsampling before the IDFT, they will fold over into your estimate of the portion of the signal that you care about. This may be undesirable. $\endgroup$ – Jason R Oct 4 '13 at 12:54
  • $\begingroup$ Jason is absolutely right. This method is a poor choice unless the data is very small in the rest of the signal compared. You will definitely get aliasing. $\endgroup$ – nispio Oct 5 '13 at 18:18
  • $\begingroup$ Yes, aliasing can be a problem, but for my problem I can track this by seeing how fast the time domain response dies down (its an impulse response) and increasing the number of points in the iFFT transform if I need longer time spans. The decay of the impulse may change from problem to problem (the iFFT is part of an optimization problem) but within one problem it will stay more or less the same. Thanks for bringing this up. $\endgroup$ – Charlie Oct 6 '13 at 2:06
5
$\begingroup$

What you're looking for is called a pruned DFT. In principle, it is possible to calculate a subset of outputs from a DFT using fewer mathematical operations. In practice, however, existing highly-optimized FFT implementations like FFTW are designed for full-output transforms. You'll find in many cases, unless you're only concerned with a very small proportion of the transform outputs, it is fastest to just calculate the entire thing and throw away what you don't want.

As it turns out, FFTW has a page in their documentation that talks about this specifically.

If you want to put in a lot of work, it is possible that you could obtain some performance gain by rolling your own FFT implementation for your specific case. FFTW is a general-purpose, any-size-transform library that works well in a variety of conditions, but it's possible that you could beat it if you tune something specifically for your application. Before expending this effort, I would think long and hard about whether it is truly a limiting factor in your system's performance.

Improve this answer
$\endgroup$
  • $\begingroup$ I did a little more digging into this, and some checking into the Pruned FFT. Wow, it's shocking. For instance in this post (mathworks.com/matlabcentral/newsreader/view_thread/91455) the OP needs to get the first 200 points only, out of 1 million points. Even with the pruned FFT the total time could be reduced a maximum of 50% but that was without using some compiler optimizations... So it looks like for my case, where I need 470 out of 16k, the full iFFT is the way to go, plus I have the extra data if I want it. Thanks for your post, I consider this "solved". $\endgroup$ – Charlie Oct 3 '13 at 19:17
  • $\begingroup$ @Charlie: If the specific set of subsamples that you are after is the first $T$ samples, then you can do better than the pruned DFT. See my answer below. $\endgroup$ – nispio Oct 3 '13 at 19:52
0
$\begingroup$

You might also look at the sparse FFT (sFFT), which, given an integer k, finds the k dominant modes in the transform of length N. It is efficient if k << N. If your data fits this problem, it might be a way to go if the modes that you are looking for are the dominant modes. It has been mentioned before: What is the sparse Fourier transform? and has its own website now: http://groups.csail.mit.edu/netmit/sFFT/index.html

If your data is complex and you are looking for a convolution, you might also look at http://fftwpp.sourceforge.net/, which both wraps FFTW into C++ for you and has implicitly dealiased convolution routines, which basically perform a pruned FFT with the appropriate padding to eliminate aliases. I'm a developer for FFTW++; if you are interested in cross-correlations on real-valued data, let us know, as there are large performance advantages to be had using implicit padding on multidemsional data.

Improve this answer
$\endgroup$
  • $\begingroup$ Could you expand on that last point? I am interested in efficient ways to find the correlation between two real signals at a small number of consecutive points. The number of points is far less than the size of the signal. $\endgroup$ – lp251 Oct 5 '13 at 19:38
  • $\begingroup$ Computing convolutions and correlations using an FFT requires zero-padding in each dimension. With implicit padding, we only expand fully in one dimension, and just re-use buffers in the other dimensions. In two or more dimensions, this saves a lot of memory, and, as it turns out, makes the algorithm much faster. Implicit padding uses a pruned FFT; we've written this for complex input with and without Hermitian symmetry, but we don't have the real-input case yet. The original paper is linked from the FFTW++ home page. $\endgroup$ – Malcolm Oct 5 '13 at 21:16
  • $\begingroup$ I looked into FFTW++ and it looked very interesting and easy to code. Unfortunately I have decided to write the code in C/C++ and compile using gcc under Cygwin (will later port to MinGW) and OpenMP doesn't seem to be supported. I will use FFTW routines and that should work for my purposes. $\endgroup$ – Charlie Oct 6 '13 at 2:03
  • $\begingroup$ FFTW++ can be compiled without OpenMP support as well; one must simply define FFTWPP_SINGLE_THREAD and it should pre-process out the OpenMP stuff. We'll definitely improve the documentation with the next release! :) $\endgroup$ – Malcolm Oct 6 '13 at 4:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.