Reply by robert bristow-johnson January 26, 20102010-01-26
On Jan 26, 6:15&#4294967295;pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 26, 9:57 am, robert bristow-johnson <r...@audioimagination.com> > wrote: > > > On Jan 25, 9:54 pm, Greg Heath <he...@alumni.brown.edu> wrote: > > > > On Jan 24, 10:30 am, robert bristow-johnson<r...@audioimagination.com> wrote: > > > > > but, even if you were doing this in straight C, with FFTW or > > > > something, you would still have to swap the halves. > > > > Because C does not allow negative indices? > > > i just realized the nature of the question. &#4294967295;the answer is "no". &#4294967295;it's > > because when you apply the window, the center of the window is where > > the phase reference should be and if you don't swap halves going into > > a regular vanilla-flavored FFT, then the center of the window is at x > > [N/2] > > I assumed that if negative indices were available swapping > would not be necessary because both fft and hanning > would have options to set the phase reference and window > center at index 0 in x(-(N-1)/2:(N-1)/2) for N odd. > > So, ... that is not possible?
sure, i s'pose. what is needed (and you could do this with MATLAB if they extended the language to include arrays having origin different than 1) would be attached attributes that would tell the calling function what the limits where (or the origin and length) for every dimension in the array. so if you had an 1xN array that had -N/2 as the index for the "first" element, then the called FFT would know it's length N and the indices run from -N/2 to -N/2 + (N-1). the no swapping would be necessary (or desired). it would be computing N/2 -1 X[k] = SUM{ x[n] * exp(-j*2*pi*n*k/N) } n=-N/2 and, unswapped, the data going in would be right. x[0] would be right in the middle and it would be treated as x[0]. but with MATLAB the way it is, if you grab a bunch of contiguous data and window it, the middle sample is x[N/2] and MATLAB treats it as such. fftshift() swaps the two halves (i was telling Jerry they should call it "swaphalves" instead of "fftshift") so that the middle sample (what used to be at x[N/2]) is treated as x[0]. but if this was done in C (or C++), an array like this would be a structure (or object) in which the length of the array and the array origin are embedded into the structure along with the contents of the array. MATLAB does half of this, which is why fft() does not need to know the length of the input array, it can infer that value. hey guys, i'm a Jerry's house right now, using Jerry's computer. this guy's a consumate engineer. he has more tools and widgets in his reality than i ever had in my imagination. r b-j
Reply by Greg Heath January 26, 20102010-01-26
On Jan 26, 9:57 am, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Jan 25, 9:54 pm, Greg Heath <he...@alumni.brown.edu> wrote: > > > On Jan 24, 10:30 am, robert bristow-johnson<r...@audioimagination.com> wrote: > > > > but, even if you were doing this in straight C, with FFTW or > > > something, you would still have to swap the halves. > > > Because C does not allow negative indices? > > i just realized the nature of the question. the answer is "no". it's > because when you apply the window, the center of the window is where > the phase reference should be and if you don't swap halves going into > a regular vanilla-flavored FFT, then the center of the window is at x > [N/2]
I assumed that if negative indices were available swapping would not be necessary because both fft and hanning would have options to set the phase reference and window center at index 0 in x(-(N-1)/2:(N-1)/2) for N odd. So, ... that is not possible? Thanks. Greg
Reply by robert bristow-johnson January 26, 20102010-01-26
On Jan 25, 9:54&#4294967295;pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 24, 10:30&#4294967295;am, robert bristow-johnson<r...@audioimagination.com> wrote: > > > but, even if you were doing this in straight C, with FFTW or > > something, you would still have to swap the halves. > > Because C does not allow negative indices?
i just realized the nature of the question. the answer is "no". it's because when you apply the window, the center of the window is where the phase reference should be and if you don't swap halves going into a regular vanilla-flavored FFT, then the center of the window is at x [N/2] which puts a phase shift of pi on all the odd samples of X[k] compared to when you *do* swap halves. you want the half with the earlier "time" indices to be at the "negative" time indices of the FFT, which is in the 2nd half (because of the inherent periodicity). if you have MATLAB cooking, try it both ways and see for yourself what the issue is, Greg. r b-j
Reply by robert bristow-johnson January 26, 20102010-01-26
On Jan 26, 6:45&#4294967295;am, Richard Dobson <richarddob...@blueyonder.co.uk>
wrote:
> On 26/01/2010 02:54, Greg Heath wrote: > .. > > >> but, even if you were doing this in straight C, with FFTW or > >> something, you would still have to swap the halves. > > > Because C does not allow negative indices? > > It does. Used all the time. Just has to result in a legal memory access.
but you have to do a teeny bit of rigermerole: double my_array[4097], *array; array = &(my_array[2048]); for (i=-2048; i<=2048; i++) { array[i] = something(); } i used to do this instead, but, of course you cannot trust that the code would always work: double neg_array[2048], array[1], pos_array[2048]; for (i=-2048; i<=2048; i++) { array[i] = something(); } if there is any justice in life, one would think that the compiler would place the three arrays in either ascending or descending order in memory and adjacent to each other. r b-j
Reply by Richard Dobson January 26, 20102010-01-26
On 26/01/2010 02:54, Greg Heath wrote:
..
>> but, even if you were doing this in straight C, with FFTW or >> something, you would still have to swap the halves. > > Because C does not allow negative indices? >
It does. Used all the time. Just has to result in a legal memory access. Richard Dobson
Reply by Greg Heath January 25, 20102010-01-25
On Jan 24, 10:30&#4294967295;am, robert bristow-johnson
<r...@audioimagination.com> wrote:
-----SNIP
> and, if MATLAB allowed us to simply define the center of the waveform > as x[0], meaning that the first half would have negative indices, and > if their FFT was smart enough to notice and roll with that, we > wouldn't be needing to to swap the two halves. > > but, even if you were doing this in straight C, with FFTW or > something, you would still have to swap the halves.
Because C does not allow negative indices?
> > 2. The centered windows squelch the data away from > > &#4294967295; &#4294967295;the center. > > dunno what that means. &#4294967295;i remember what "squelch" means from my old > ham days (but it was the CB radios that had such a knob).
I should have used :attenuate" .
> > &#4294967295; &#4294967295;Therefore the window shouldn't > > &#4294967295; &#4294967295;slide more than half of its length. > > with a symmetric window, it should slide exactly half the length. > > > I am still bothered by the coherent addition of data > > with different phase references; i.e.should it be > > coherent addition or noncoherent averaging? > > well, maybe in the frequency domain, one needs to be careful to keep > things phase coherent. &#4294967295;this is what we do in a phase vocoder, for > example.
OK. Thanks. Greg
Reply by robert bristow-johnson January 24, 20102010-01-24
On Jan 24, 10:30&#4294967295;am, robert bristow-johnson
<r...@audioimagination.com> wrote:
> On Jan 24, 6:05&#4294967295;am, Greg Heath <he...@alumni.brown.edu> wrote:
...
> > &#4294967295; &#4294967295;Therefore the window shouldn't > > &#4294967295; &#4294967295;slide more than half of its length. > > with a symmetric window, it should slide exactly half the length. > > > I am still bothered by the coherent addition of data > > with different phase references; i.e.should it be > > coherent addition or noncoherent averaging? > > well, maybe in the frequency domain, one needs to be careful to keep > things phase coherent. &#4294967295;this is what we do in a phase vocoder, for > example.
i see that MATLAB changed the name from the incorrect "hanning" to the correct "hann", but, if i recall, they might be a little funky about the definition. so setting aside a possible inexactness, i mean to (at least for even N) set it up so that the adjacent hann windows were perfectly complimentary in how they added. so, at least for the null processing (Y=X), it should be clear how things remain phase coherent. if things are not so coherent for some non-null processing, then maybe someone's frequency-domain processor could use a little more work. r b-j
Reply by robert bristow-johnson January 24, 20102010-01-24
On Jan 24, 6:05&#4294967295;am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 24, 4:44 am, Greg Heath <he...@alumni.brown.edu> wrote: > > > > > On Jan 23, 1:20 pm, robert bristow-johnson <r...@audioimagination.com> > > wrote: > -----SNIP > > > so does this make sense?: > > > > &#4294967295; &#4294967295; y_output = zeros(size(x_input)); > > > &#4294967295; &#4294967295; N = length(x); > > > &#4294967295; &#4294967295; framehop = floor(N/2); > > > &#4294967295; &#4294967295; n = 0; > > > &#4294967295; &#4294967295; while (n < length(x_input)-N) > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; x = x_input(n+1:n+N); > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; x = x .* hanning(size(x)); > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; X = fft(ifftshift(x)); > > > % > > > % &#4294967295; &#4294967295; &#4294967295; do your processing of X here. > > > % &#4294967295; &#4294967295; &#4294967295; phase of X is nicely ironed out. > > > % > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; Y = X; &#4294967295; &#4294967295; &#4294967295; % null processing > > > % > > > % > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; y = fftshift(ifft(Y)); > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; y_output(n+1:n+N) = y_output(n+1:n+N) + y; > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; n = n + framehop; > > > &#4294967295; &#4294967295; end; > > > > would you say that's correct, Greg, or should the first > > > one be fftshift and the second be ifftshift? > > > The order of ifftshift=>fft=>ifft=>fftshift is correct. > > Therefore disproving my recommendation of only using > > iffshift and fftshift in the temporal and spectral > > domains, respectively. > > > However, although the phase of X is referenced to the > > middle of x, DC is still at the left edge of X. Chances > > are the processing referred to involves the assumption > > that the spectrum is zero centered. For that scenario > >change the code as follows > > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; X = fft(ifftshift(x)); > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;Xb = fftshift(X) > > > % &#4294967295; &#4294967295; &#4294967295; do your processing of Xb here. > > % &#4294967295; &#4294967295; &#4294967295; phase of Xb is nicely ironed out. > > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;Yb = function(Xb); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; % processing > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; Y &#4294967295;= ifftshift(Yb); > > > On the other hand, I do not comprehend what the algorithm > > is doing. You are sliding the window 1/2 of it's length and > > coherently adding the results of the next iteration. > > > Should the addition be noncoherent vis a vis Welch? > > Otherwise I don't understand.. > > > > and do you see why one might want to use one or the other > > > on the input to the fft? > > > Not yet. I'm still staring at the headlights. > > > > without considering if further, i am not so sure how the overlap- > > > adding will work for odd N. &#4294967295;but it should work transparently (except > > > for the two ends of y_output) for even N. &#4294967295;maybe it should be framehop > > > = ceil(N/2), but that doesn't seem right either. > > > Coherent overlap adding is a new one on me. What is the purpose? >
phase vocoder. or any generic kind of frequency-domain processing (but don't confuse with the overlap-add or overlap-save fast FIR filters, there is some difference).
> OK. I see part of it now. > > 1. The MATLAB windows are centered at the middle of > &#4294967295; &#4294967295;waveforms.
whether it's MATLAB or not.
> Therefore, that is where you want > &#4294967295; &#4294967295;the phase reference.
and, if MATLAB allowed us to simply define the center of the waveform as x[0], meaning that the first half would have negative indices, and if their FFT was smart enough to notice and roll with that, we wouldn't be needing to to swap the two halves. but, even if you were doing this in straight C, with FFTW or something, you would still have to swap the halves.
> 2. The centered windows squelch the data away from > &#4294967295; &#4294967295;the center.
dunno what that means. i remember what "squelch" means from my old ham days (but it was the CB radios that had such a knob).
> Therefore the window shouldn't > &#4294967295; &#4294967295;slide more than half of its length.
with a symmetric window, it should slide exactly half the length.
> I am still bothered by the coherent addition of data > with different phase references; i.e.should it be > coherent addition or noncoherent averaging?
well, maybe in the frequency domain, one needs to be careful to keep things phase coherent. this is what we do in a phase vocoder, for example. r b-j
Reply by Greg Heath January 24, 20102010-01-24
On Jan 24, 4:44 am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 23, 1:20 pm, robert bristow-johnson <r...@audioimagination.com> > wrote:
-----SNIP
> > so does this make sense?: > > > y_output = zeros(size(x_input)); > > N = length(x); > > framehop = floor(N/2); > > n = 0; > > while (n < length(x_input)-N) > > x = x_input(n+1:n+N); > > x = x .* hanning(size(x)); > > X = fft(ifftshift(x)); > > % > > % do your processing of X here. > > % phase of X is nicely ironed out. > > % > > Y = X; % null processing > > % > > % > > y = fftshift(ifft(Y)); > > y_output(n+1:n+N) = y_output(n+1:n+N) + y; > > n = n + framehop; > > end; > > > would you say that's correct, Greg, or should the first > > one be fftshift and the second be ifftshift? > > The order of ifftshift=>fft=>ifft=>fftshift is correct. > Therefore disproving my recommendation of only using > iffshift and fftshift in the temporal and spectral > domains, respectively. > > However, although the phase of X is referenced to the > middle of x, DC is still at the left edge of X. Chances > are the processing referred to involves the assumption > that the spectrum is zero centered. For that scenario >change the code as follows > > > X = fft(ifftshift(x)); > > Xb = fftshift(X) > > % do your processing of Xb here. > % phase of Xb is nicely ironed out. > > Yb = function(Xb); % processing > Y = ifftshift(Yb); > > On the other hand, I do not comprehend what the algorithm > is doing. You are sliding the window 1/2 of it's length and > coherently adding the results of the next iteration. > > Should the addition be noncoherent vis a vis Welch? > Otherwise I don't understand.. > > > and do you see why one might want to use one or the other > > on the input to the fft? > > Not yet. I'm still staring at the headlights. > > > without considering if further, i am not so sure how the overlap- > > adding will work for odd N. but it should work transparently (except > > for the two ends of y_output) for even N. maybe it should be framehop > > = ceil(N/2), but that doesn't seem right either. > > Coherent overlap adding is a new one on me. What is the purpose?
OK. I see part of it now. 1. The MATLAB windows are centered at the middle of waveforms. Therefore, that is where you want the phase reference. 2. The centered windows squelch the data away from the center. Therefore the window shouldn't slide more than half of its length. I am still bothered by the coherent addition of data with different phase references; i.e.should it be coherent addition or noncoherent averaging? Hope this helps. Greg
Reply by Greg Heath January 24, 20102010-01-24
On Jan 23, 1:20 pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Jan 23, 4:55 am, Greg Heath <he...@alumni.brown.edu> wrote: > > On Jan 21, 10:26 pm, robert bristow-johnson > > > <r...@audioimagination.com> wrote: > > > On Jan 21, 3:48 pm, Jerry Avins <j...@ieee.org> wrote: > > > > robert bristow-johnson wrote: > > > > > so, if you're gonna build a case that ifftshift() should be used only > > > > > with one and that fftshift() should be used only with the other, it's > > > > > already a loser. > > > > > Greg's claim is that for naive users who need one simple rule that > > > > always works properly, it's a good way to go. That seems plausible. > > > > i just dunno if the difference between fftshift and ifftshift has > > > anything to do with the difference between fft and ifft (of which the > > > latter has no qualitative difference other than scaling). > > > An important difference is that f = 0 has physical significance > > whereas t = 0 does not. > > that's making an assumption of the physical meaning of x[] or X[]. i > consider the DFT and iDFT to be qualitatively equivalant. > > > No matter how x(t) is circularly shifted, > > X = fft(x) always yields X(1) = X(f=0). > > you can say the same about the inverse. no matter how X(f) is > circularly shifted, x = ifft(X) always yields x(1) = x(t=0). > > > I'm not familiar with the source code of fftshift for > > other languages. > > i've never seen a library with a function named "fftshift" other than > MATLAB. > > > > in the "time domain", i dunno whether it's better, for mathematical > > > reasons, to map x[(N+1)/2] to x[0] or to map X[(N-1)/2] to x[0]. i > > > think maybe it's the latter. so i think it's not better to use > > > ifftshift to fix the phasey issue than fftshift. but i am not > > > entirely sure. > > > The important point is, for symmetric time functions, to map the > > symmetry point to x(1). > > > fftshift([ 2 1 0 1 2 ]) = [1 2 2 1 0] > > > wheras > > > fftshift([ 2 1 0 1 2 ]) = [0 1 2 2 1] > > there's a missing "i",
WHOOPS! That should be "i"fftshift.
> so let's see what you meant. quoting from before: > > ifftshift(x(1:N)) = [ x(1+floor(N/2):N), x(1:floor(N/2)) ] > fftshift(X(1:N)) = [ X(1+ceil(N/2):N), X(1:ceil(N/2)) ] > > i think it's > > fftshift([ 2 1 0 1 2 ]) = [1 2 2 1 0] > > wheras > > ifftshift([ 2 1 0 1 2 ]) = [0 1 2 2 1] > > is that not correct?
Yep.
> > Since it is desirable that the same function be used > > when N is even or odd and large or small, the > > choice is pretty much limited to the approach I > > have recommended. > > for the purpose of mapping the point of symmetry to x(1) (what the > rest of the world calls x[0]), then you should use ifftshift for odd N > (or even N).
Yes.
> so does this make sense?: > > y_output = zeros(size(x_input)); > N = length(x); > framehop = floor(N/2); > n = 0; > while (n < length(x_input)-N) > x = x_input(n+1:n+N); > x = x .* hanning(size(x)); > X = fft(ifftshift(x)); > % > % do your processing of X here. > % phase of X is nicely ironed out. > % > Y = X; % null processing > % > % > y = fftshift(ifft(Y)); > y_output(n+1:n+N) = y_output(n+1:n+N) + y; > n = n + framehop; > end; > > would you say that's correct, Greg, or should the first one be > fftshift and the second be ifftshift?
The order of ifftshift=>fft=>ifft=>fftshift is correct. Therefore disproving my recommendation of only using iffshift and fftshift in the temporal and spectral domains, respectively. However, although the phase of X is referenced to the middle of x, DC is still at the left edge of X. Chances are the processing referred to involves the assumption that the spectrum is zero centered. For that scenario change the code as follows
> X = fft(ifftshift(x));
Xb = fftshift(X) % do your processing of Xb here. % phase of Xb is nicely ironed out. Yb = function(Xb); % processing Y = ifftshift(Yb); On the otherhand, I do not comprehend what the algorithm is doing. You are sliding the window 1/2 of it's length and coherently adding the results of the next iteration. Should the addition be noncoherent vis a vis Welch? Otherwise I don't understand..
> and do you see why one might want to use one or the other > on the input to the fft?
Not yet. I'm still staring at the headlights.
> without considering if further, i am not so sure how the overlap- > adding will work for odd N. but it should work transparently (except > for the two ends of y_output) for even N. maybe it should be framehop > = ceil(N/2), but that doesn't seem right either.
Coherent overlap adding is a new one on me. What is the purpose? Hope this helps. Greg