Forums

FFTW convolution scaling problem

Started by schoopy January 22, 2008
Hello everyone,

I need help on this one.
I am trying to use FFTW to compute the integral of x^-4 and some smooth
on/off signal first in 1D. But already here I am experiencing some scaling
problem by a factor close to 3.
This is what I do (after reading bunch of web pages) for a N sampling:
1. compute N points of the functions between -L and L.
2. wrap around the response function (x^-4) in a 2N vector as such:
positive x, zero, negative x values
3. zero pad the signal by placing it between N/4 and 3N/4.
4. do the FFTs
5. scale by a factor 1/(2N)
6. complex multiply the results.
7. do the iFFT
8. scale by 2L.

The shape of my result is perfect except a residual scale factor around 3
(2.96...). 
So far I tested changing N or L but results remain stable. I cannot
explain this but it seems someone add similar problem, the answer is not
published though!

Any idea is welcome !!!

Thanks


Have you resolve scaling problem? I got the same issue.


On 6/29/12 9:03 PM, azz750 wrote:
> Have you resolve scaling problem? I got the same issue. >
okay, azz, this is the second time you've posted about this. perhaps nobody responded (until now) because you haven't posed or formulated the question with sufficient description that we have any idea what you're typing about. i have an idea about scaling and the DFT (which is what FFTW does). but i don't know of a "scaling problem" regarding FFTW. is it the lack of a 1/N in the inverse transform? -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On 6/30/2012 12:35 AM, robert bristow-johnson wrote:
> On 6/29/12 9:03 PM, azz750 wrote: >> Have you resolve scaling problem? I got the same issue. >> >
> > okay, azz, this is the second time you've posted about this. perhaps > nobody responded (until now) because you haven't posed or formulated the > question with sufficient description that we have any idea what you're > typing about. > > i have an idea about scaling and the DFT (which is what FFTW does). but > i don't know of a "scaling problem" regarding FFTW. is it the lack of a > 1/N in the inverse transform? > >
Good point. May be OP is talking about similar question for Matlab's fft()? Where on the help page, there is an example that divides the fft() result by length of the signal. Also it multiplies the abs(fft) by 2 later on. These are scaling factors. I've seen few questions on these scaling factors at the matlab news group. Here is the Mathworks web page with this example http://www.mathworks.com/help/techdoc/ref/fft.html if you look at the example shown there, you'll see these 2 lines Y = fft(y,NFFT)/L; % notice the divide and plot(f,2*abs(Y(1:NFFT/2+1))) % notice the 2 May be the OP is asking about these? just guessing. ps. These have been talked about and answered many times in the Matlab newsgroup. May be OP can google this if this is what they are asking about. --Nasser
On 6/29/2012 8:03 PM, azz750 wrote:
> Have you resolve scaling problem? I got the same issue. > >
If I recall correctly, the FFTW functions do not scale the output, neither the forward or the inverse. Any scaling is up to you. OUP
On 6/30/12 9:00 PM, One Usenet Poster wrote:
> On 6/29/2012 8:03 PM, azz750 wrote: >> Have you resolve scaling problem? I got the same issue. >> >> > > If I recall correctly, the FFTW functions do not scale the output, > neither the forward or the inverse. Any scaling is up to you.
that's what i thought, and i tried to find the document that said so, and i didn't. but i think that's the case. if you do FFT (W) followed by iFFT, you will get your original data scaled up by a factor of N. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Indeed me and topic starter know about scaling by N which is required by fftw to get original function after forward and inverse transformations (hence "5. scale by a factor 1/(2N)" in original question).

I had following problem:

- Calculate iFFT ( FFT(function)^n ), here function is defined in d-dimensional cube with side L, n = 1, 2, ....

I tried to perform this calculation with Gaussian as a "function" in d=2. This is because it can be done analytically (btw, the answer is the Gauss with wider "width" and some factor in front). However, calculations done with fftw gave me very different factor (although, the width was the same).

Indeed the problem is "8. scale by 2L" I found that in my case the result of calculations done with fftw should be scaled by ((2L/N)^d)^(n-1)/(N^d)

--http://compgroups.net/comp.dsp/fftw-convolution-scaling-problem/1176930


On Jun 30, 8:00=A0pm, One Usenet Poster <m...@my.computer.org> wrote:
> On 6/29/2012 8:03 PM, azz750 wrote: > > > Have you resolve scaling problem? I got the same issue. > > If I recall correctly, the FFTW functions do not scale the output, > neither the forward or the inverse. =A0Any scaling is up to you.
I haven't studied FFTW in depth (just a brief review), so I can't say too much about what it's doing. But I discussed the use of the complex fourier transform along with its real forms (in particular the cosine transform) for use with convolution and posted some demo code for convolution-by-complex fourier transform. Notes on the DFT, Convolution and Implementations http://www.docstoc.com/docs/124401006/Notes-on-the-DFT-Convolution-and-Impl= ementationsThe inverse should have the divide by N. The basic routine for use with the complex fourier transform is: Convolve sequences A and B of lengths respectively M and N: (1) Find a suitable size P (e.g. power of 2) at least as large as M + N - 1, where (2) Zero-pad A on the last P - M slots and transform it to A' (3) Zero-pad B on the last P - N slots and transform it to B' (4) Multiply A' and B' point-wise to get C' (also of length P) (5) Inverse transform C' to get C. Then truncate to length M + N - 1. The routine for use with the cosine transform: Convolve sequences A and B of lengths respectively M and N: (0) If M is even, N odd, switch the two sequences (1) Find a suitable number P at least as large as 2(M + N) - (N + 1) mod 2. (2) Zero-pad A on the LEFT by (N + 1) div 2, and on the right to length P. Transform to A'. (3) Zero-pad B on the LEFT by (N + 1) div 2 + M, and on the right to length P. Transform to B'. (4) Multiply A', B' point-wise to get C'. (5) Inverse transform C' to C, shift to the left by M + N + (N mod 2) and truncate on the right to length M + N - 1. The explicit forms for the transforms are given in the PDF (along with the scaling factors used for the inverses). The reasons for all the shiftings are explained in the PDF. The same general analysis can be applied to the other real-form transforms, though I didn't carry out any of the analyses.