Re: [Numpy-discussion] Bug or surprising undocumented behaviour in irfft
On 8/29/07, Anne Archibald <[EMAIL PROTECTED]> wrote: > > On 29/08/2007, Charles R Harris <[EMAIL PROTECTED]> wrote: > > > > Is this also appropriate for the other FFTs? (inverse real, complex, > > > hermitian, what have you) I have written a quick hack (attached) that > > > should do just that rescaling, but I don't know that it's a good idea, > > > as implemented. Really, for a complex IFFT it's extremely peculiar to > > > add the padding where we do (between frequency -1 and frequency zero); > > > it would make more sense to pad at the high frequencies (which are in > > > the middle of the array). Forward FFTs, though, can reasonably be > > > padded at the end, and it doesn't make much sense to rescale the last > > > data point. > > > > It all depends on the data and what you intend. Much of my experience is > > with Michaelson interferometers and in that case the interferogram is > > essentially an autocorrelation, so it is desirable to keep its center at > > sample zero and let the left side wrap around, so ideally you fill in > the > > middle as you suggest. You can also pad at the end if you don't put the > > center at zero, but then you need to phase shift the spectrum in a way > that > > corresponds to rotating the center to index zero and padding in the > middle. > > I expect you would want to do the same thing for complex transforms if > they > > are of real data and do the nyquist divided by two thingy. If the high > > frequencies in a complex transform are actually high frequencies and not > > aliases of negative frequencies, then you will want to just append > zeros. > > That case also occurs, I have designed decimating complex filters that > > produce output like that, they were like single sideband in the radi o > > world. > > So is it a fair summary to say that for irfft, it is fairly clear that > one should adjust the Nyquist coefficient, but for the other varieties > of FFT, the padding done by numpy is just one of many possible > choices? > > Should numpy be modified so that irfft adjusts the Nyquist > coefficient? Should this happen only for irfft? Yes, I think that should be the case. If the complex transforms pad in the middle, then they are treating the high frequencies as aliases, but unless they explicitly duplicate the Nyquist coefficient scaling isn't needed. Hmm, actually, I think that is wrong. The original data points will be reproduced, but what happens in between points? In between there is a difference between positive and negative frequences. So in a complex transform of real data one would want to split the Nyquist coefficient between high and low frequencies. I don't think it is possible to make a general statement about the complex case. Just hope the middle frequency is zero so you can ignore the problem ;) What happens in the real case is that the irfft algorithm uses the Hermitean symmetry of the spectrum, so the coefficient is implicitly duplicated. > I usually multiply the forward transform by the sample interval, in secs > or > > cm, and the unscaled inverse transform by the frequency sample interval, > in > > Hz or cm^-1. That treats both the forward and inverse fft like > > approximations to the integral transforms and makes the units those of > > spectral density. If you think trapezoidal rule, then you will also see > > factors of .5 at the ends, but that is a sort of apodization that is > > consistent with how Fourier series converge at discontinuities. In the > > normal case where no interpolation is done the product of the sample > > intervals is 1/N, so it reduces to the usual convention. Note that in > your > > example the sampling interval decreases when you do the interpolation, > so if > > you did another forward transform it would be scaled down to account for > the > > extra points in the data. > > That's a convenient normalization. > > Do you know if there's a current package to associate units with numpy > arrays? For my purposes it would usually be sufficient to have arrays > of quantities with uniform units. Conversions need only be > multiplicative (I don't care about Celsius-to-Fahrenheit style > conversions) and need not even be automatic, though of course that > would be convenient. Right now I use Frink for that sort of thing, but > it would have saved me from making a number of minor mistakes in > several pieces of python code I've written. There was a presentation by some fellow from CalTech at SciPy 2005 (4?) about such a system, but ISTR it looked pretty complex. C++ template programming does it with traits and maybe the Enthought folks have something useful along those lines. Otherwise, I don't know of any such system for general use. Maybe ndarray could be subclassed? It can be convenient to multiply and divide units, so maybe some sort of string with something to gather the same units together with a power could be a useful way to track them and wouldn't tie one down to any particular choice. Chuck ___
Re: [Numpy-discussion] Bug or surprising undocumented behaviour in irfft
On 29/08/2007, Charles R Harris <[EMAIL PROTECTED]> wrote: > > Is this also appropriate for the other FFTs? (inverse real, complex, > > hermitian, what have you) I have written a quick hack (attached) that > > should do just that rescaling, but I don't know that it's a good idea, > > as implemented. Really, for a complex IFFT it's extremely peculiar to > > add the padding where we do (between frequency -1 and frequency zero); > > it would make more sense to pad at the high frequencies (which are in > > the middle of the array). Forward FFTs, though, can reasonably be > > padded at the end, and it doesn't make much sense to rescale the last > > data point. > > It all depends on the data and what you intend. Much of my experience is > with Michaelson interferometers and in that case the interferogram is > essentially an autocorrelation, so it is desirable to keep its center at > sample zero and let the left side wrap around, so ideally you fill in the > middle as you suggest. You can also pad at the end if you don't put the > center at zero, but then you need to phase shift the spectrum in a way that > corresponds to rotating the center to index zero and padding in the middle. > I expect you would want to do the same thing for complex transforms if they > are of real data and do the nyquist divided by two thingy. If the high > frequencies in a complex transform are actually high frequencies and not > aliases of negative frequencies, then you will want to just append zeros. > That case also occurs, I have designed decimating complex filters that > produce output like that, they were like single sideband in the radi o > world. So is it a fair summary to say that for irfft, it is fairly clear that one should adjust the Nyquist coefficient, but for the other varieties of FFT, the padding done by numpy is just one of many possible choices? Should numpy be modified so that irfft adjusts the Nyquist coefficient? Should this happen only for irfft? > I usually multiply the forward transform by the sample interval, in secs or > cm, and the unscaled inverse transform by the frequency sample interval, in > Hz or cm^-1. That treats both the forward and inverse fft like > approximations to the integral transforms and makes the units those of > spectral density. If you think trapezoidal rule, then you will also see > factors of .5 at the ends, but that is a sort of apodization that is > consistent with how Fourier series converge at discontinuities. In the > normal case where no interpolation is done the product of the sample > intervals is 1/N, so it reduces to the usual convention. Note that in your > example the sampling interval decreases when you do the interpolation, so if > you did another forward transform it would be scaled down to account for the > extra points in the data. That's a convenient normalization. Do you know if there's a current package to associate units with numpy arrays? For my purposes it would usually be sufficient to have arrays of quantities with uniform units. Conversions need only be multiplicative (I don't care about Celsius-to-Fahrenheit style conversions) and need not even be automatic, though of course that would be convenient. Right now I use Frink for that sort of thing, but it would have saved me from making a number of minor mistakes in several pieces of python code I've written. Thanks, Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Bug or surprising undocumented behaviour in irfft
Hi Anne, On 8/29/07, Anne Archibald <[EMAIL PROTECTED]> wrote: > > On 29/08/2007, Charles R Harris <[EMAIL PROTECTED]> wrote: > > > > > What is going on is that the coefficient at the Nyquist frequency > appears > > once in the unextended array, but twice when the array is extended with > > zeros because of the Hermitean symmetry. That should probably be fixed > in > > the upsampling code. > > Is this also appropriate for the other FFTs? (inverse real, complex, > hermitian, what have you) I have written a quick hack (attached) that > should do just that rescaling, but I don't know that it's a good idea, > as implemented. Really, for a complex IFFT it's extremely peculiar to > add the padding where we do (between frequency -1 and frequency zero); > it would make more sense to pad at the high frequencies (which are in > the middle of the array). Forward FFTs, though, can reasonably be > padded at the end, and it doesn't make much sense to rescale the last > data point. It all depends on the data and what you intend. Much of my experience is with Michaelson interferometers and in that case the interferogram is essentially an autocorrelation, so it is desirable to keep its center at sample zero and let the left side wrap around, so ideally you fill in the middle as you suggest. You can also pad at the end if you don't put the center at zero, but then you need to phase shift the spectrum in a way that corresponds to rotating the center to index zero and padding in the middle. I expect you would want to do the same thing for complex transforms if they are of real data and do the nyquist divided by two thingy. If the high frequencies in a complex transform are actually high frequencies and not aliases of negative frequencies, then you will want to just append zeros. That case also occurs, I have designed decimating complex filters that produce output like that, they were like single sideband in the radio world. > The inverse irfft also scales by dividing by the new transform size > instead > > of the original size, so the result needs to be scaled up for the > > interpolation to work. It is easy to go wrong with fft's when the > correct > > sampling/frequency scales aren't carried with the data. I always do that > > myself so that the results are independent of transform > size/interpolation > > and expressed in some standard units. > > The scaling of the FFT is a pain everywhere. I always just try it a > few times until I get the coefficients right. I sort of like FFTW's > convention of never normalizing anything - it means the transforms > have nice simple formulas, though unfortunately it also means that > ifft(fft(A))!=A. In any case the normalization of numpy's FFTs is not > something that can reasonably be changed, even in the special case of > the zero-padding inverse (and forward) FFTs. I usually multiply the forward transform by the sample interval, in secs or cm, and the unscaled inverse transform by the frequency sample interval, in Hz or cm^-1. That treats both the forward and inverse fft like approximations to the integral transforms and makes the units those of spectral density. If you think trapezoidal rule, then you will also see factors of .5 at the ends, but that is a sort of apodization that is consistent with how Fourier series converge at discontinuities. In the normal case where no interpolation is done the product of the sample intervals is 1/N, so it reduces to the usual convention. Note that in your example the sampling interval decreases when you do the interpolation, so if you did another forward transform it would be scaled down to account for the extra points in the data. Chuck ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Bug or surprising undocumented behaviour in irfft
On 29/08/2007, Charles R Harris <[EMAIL PROTECTED]> wrote: > > > What is going on is that the coefficient at the Nyquist frequency appears > once in the unextended array, but twice when the array is extended with > zeros because of the Hermitean symmetry. That should probably be fixed in > the upsampling code. Is this also appropriate for the other FFTs? (inverse real, complex, hermitian, what have you) I have written a quick hack (attached) that should do just that rescaling, but I don't know that it's a good idea, as implemented. Really, for a complex IFFT it's extremely peculiar to add the padding where we do (between frequency -1 and frequency zero); it would make more sense to pad at the high frequencies (which are in the middle of the array). Forward FFTs, though, can reasonably be padded at the end, and it doesn't make much sense to rescale the last data point. > The inverse irfft also scales by dividing by the new transform size instead > of the original size, so the result needs to be scaled up for the > interpolation to work. It is easy to go wrong with fft's when the correct > sampling/frequency scales aren't carried with the data. I always do that > myself so that the results are independent of transform size/interpolation > and expressed in some standard units. The scaling of the FFT is a pain everywhere. I always just try it a few times until I get the coefficients right. I sort of like FFTW's convention of never normalizing anything - it means the transforms have nice simple formulas, though unfortunately it also means that ifft(fft(A))!=A. In any case the normalization of numpy's FFTs is not something that can reasonably be changed, even in the special case of the zero-padding inverse (and forward) FFTs. Anne fftfix Description: Binary data ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Bug or surprising undocumented behaviour in irfft
On 8/29/07, Charles R Harris <[EMAIL PROTECTED]> wrote: > > Anne, > > On 8/29/07, Anne Archibald <[EMAIL PROTECTED]> wrote: > > > > Hi, > > > > numpy's Fourier transforms have the handy feature of being able to > > upsample and downsample signals; for example the documentation cites > > irfft(rfft(A),16*len(A)) as a way to get a Fourier interpolation of A. > > However, there is a peculiarity with the way numpy handles the > > highest-frequency coefficient. > > > > > The upshot is, if I correctly understand what is going on, that the > > last coefficient needs to be treated somewhat differently from the > > others; when one pads with zeros in order to upsample the signal, one > > should multiply the last coefficient by 0.5. Should this be done in > > numpy's upsampling code? Should it at least be documented? > > > What is going on is that the coefficient at the Nyquist frequency appears > once in the unextended array, but twice when the array is extended with > zeros because of the Hermitean symmetry. That should probably be fixed in > the upsampling code. > The inverse irfft also scales by dividing by the new transform size instead of the original size, so the result needs to be scaled up for the interpolation to work. It is easy to go wrong with fft's when the correct sampling/frequency scales aren't carried with the data. I always do that myself so that the results are independent of transform size/interpolation and expressed in some standard units. In [9]: a = array([1, 0, 0, 0], dtype=double) In [10]: b = rfft(a) In [11]: b[2] *= .5 In [12]: irfft(b,8) Out[12]: array([ 0.5 , 0.3017767, 0. , -0.0517767, 0. , -0.0517767, 0. , 0.3017767]) In [13]: 2*irfft(b,8) Out[13]: array([ 1., 0.60355339, 0., -0.10355339, 0., -0.10355339, 0., 0.60355339]) I don't know where that should be fixed. Chuck ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Bug or surprising undocumented behaviour in irfft
Anne, On 8/29/07, Anne Archibald <[EMAIL PROTECTED]> wrote: > > Hi, > > numpy's Fourier transforms have the handy feature of being able to > upsample and downsample signals; for example the documentation cites > irfft(rfft(A),16*len(A)) as a way to get a Fourier interpolation of A. > However, there is a peculiarity with the way numpy handles the > highest-frequency coefficient. The upshot is, if I correctly understand what is going on, that the > last coefficient needs to be treated somewhat differently from the > others; when one pads with zeros in order to upsample the signal, one > should multiply the last coefficient by 0.5. Should this be done in > numpy's upsampling code? Should it at least be documented? What is going on is that the coefficient at the Nyquist frequency appears once in the unextended array, but twice when the array is extended with zeros because of the Hermitean symmetry. That should probably be fixed in the upsampling code. Chuck ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Bug or surprising undocumented behaviour in irfft
Hi, numpy's Fourier transforms have the handy feature of being able to upsample and downsample signals; for example the documentation cites irfft(rfft(A),16*len(A)) as a way to get a Fourier interpolation of A. However, there is a peculiarity with the way numpy handles the highest-frequency coefficient. First of all, the normalization: In [65]: rfft(cos(2*pi*arange(8)/8.)) Out[65]: array([ -3.44505240e-16 +0.e+00j, 4.e+00 -1.34392280e-15j, 1.22460635e-16 -0.e+00j, -1.16443313e-16 -8.54080261e-16j, 9.95839695e-17 +0.e+00j]) In [66]: rfft(cos(2*4*pi*arange(8)/8.)) Out[66]: array([ 0.+0.j, 0.+0.j, 0.-0.j, 0.+0.j, 8.+0.j]) So a cosine signal gives 0.5*N if its frequency F is 0http://projects.scipy.org/mailman/listinfo/numpy-discussion
