Hi all, Half a year I posted here about an inconsistency between DFT and its continuous repeated, interpolated Fourier series which should match (up to a factor). I did not really get an answer to it. So I try it again - step by step. My main point is that this is not true except we do NOT use the interplated signal but rather the infinite-bandwidth, non-interpolated signal! So I quote the relevant sentences from https://ccrma.stanford.edu/~jos/mdft/Relation_DFT_Fourier_Series.html (or http://www.dsprelated.com/dspbooks/mdft/Relation_DFT_Fourier_Series.html) and use "JOS>> " as ident. JOS>> We now show that the DFT of a sampled signal x(n) (of length N) JOS>> ... Ok, I plotted one: x = [ 2 1 ], N=2: http://snag.gy/Q4AWO.jpg Let's call it "x[n]". JOS>> ... is proportional to the Fourier series coefficients of the JOS>> continuous periodic signal obtained by repeating ... Ok, converted the signal to a continuous signal and repeated: http://snag.gy/woEyc.jpg Let's call it "xs(t)". Note that: 1.) I zoomed into an arbitrary time window (it's periodic) 2.) The arrows represent weighted Diracs 3.) The (continuous) values between the Diracs are zero(!) 4.) The signal has infinite bandwidth 5.) The period T0 = 2 JOS>> ... and interpolating x. Ok, interpolated: http://snag.gy/TsfPr.jpg Let's call it "x(t)". For better clarity, the non-interpolated xs(t) is also included. Note that the interpolation corresponds to the convolution with an appropriately weighted sinc or equivalently, filtering with the appropriate box lowpass filter. *Now* the analog signal is bandlimited. And it is *nonzero* everywhere! And according to the quoted sentence (I read it about thousand times) the Fourier coefficients obtained from x(t) over one period should be proportional to the DFT coefficients of x[n]. Then he does a bunch of stuff, but then he writes: JOS>> We wish to find the continuous-time Fourier series of the JOS>> sampled periodic signal x(nT) NO!! You told something else above! We wish to find the "continuous-time Fourier series of the periodic, *interpolated* signal x(t)" If I continue this way - of course, the result is obvious because I multiply with the Dirac comb and I am back to xs(t). But that's not my bandlimited, analog version of x[n]. I delete the non-zero data which was introduced by the interpolation. I know that the relation is conceptionally correct but I insist that the quoted document is not written consistently (if you do not agree could you show why it is correct?) In any case, how can I find the relation between DFT{x[n]} and the Fourier series of x(t)? Thanks Peter
Relation of the DFT to Fourier Series revisited
Started by ●August 6, 2014
Reply by ●August 6, 20142014-08-06
On Wed, 06 Aug 2014 12:31:01 -0700, Peter Mairhofer wrote:> Hi all, > > Half a year I posted here about an inconsistency between DFT and its > continuous repeated, interpolated Fourier series which should match (up > to a factor). I did not really get an answer to it. So I try it again - > step by step. > > My main point is that this is not true except we do NOT use the > interplated signal but rather the infinite-bandwidth, non-interpolated > signal! > > So I quote the relevant sentences from > https://ccrma.stanford.edu/~jos/mdft/Relation_DFT_Fourier_Series.html > (or > http://www.dsprelated.com/dspbooks/mdft/Relation_DFT_Fourier_Series.html)> and use "JOS>> " as ident. > > JOS>> We now show that the DFT of a sampled signal x(n) (of length N) > JOS>> ... > > Ok, I plotted one: x = [ 2 1 ], N=2: > > http://snag.gy/Q4AWO.jpg > > Let's call it "x[n]". > > JOS>> ... is proportional to the Fourier series coefficients of the > JOS>> continuous periodic signal obtained by repeating ... > > Ok, converted the signal to a continuous signal and repeated: > > http://snag.gy/woEyc.jpg > > Let's call it "xs(t)". > > Note that: > 1.) I zoomed into an arbitrary time window (it's periodic) > 2.) The arrows represent weighted Diracs 3.) The (continuous) values > between the Diracs are zero(!) > 4.) The signal has infinite bandwidth 5.) The period T0 = 2 > > JOS>> ... and interpolating x. > > Ok, interpolated: > > http://snag.gy/TsfPr.jpg > > Let's call it "x(t)". > > For better clarity, the non-interpolated xs(t) is also included. > > Note that the interpolation corresponds to the convolution with an > appropriately weighted sinc or equivalently, filtering with the > appropriate box lowpass filter. > > *Now* the analog signal is bandlimited. > And it is *nonzero* everywhere! > > And according to the quoted sentence (I read it about thousand times) > the Fourier coefficients obtained from x(t) over one period should be > proportional to the DFT coefficients of x[n]. > > Then he does a bunch of stuff, but then he writes: > > JOS>> We wish to find the continuous-time Fourier series of the JOS>> > sampled periodic signal x(nT) > > NO!! You told something else above! We wish to find the "continuous-time > Fourier series of the periodic, *interpolated* signal x(t)" > > If I continue this way - of course, the result is obvious because I > multiply with the Dirac comb and I am back to xs(t). But that's not my > bandlimited, analog version of x[n]. I delete the non-zero data which > was introduced by the interpolation. > > I know that the relation is conceptionally correct but I insist that the > quoted document is not written consistently (if you do not agree could > you show why it is correct?) > > In any case, how can I find the relation between DFT{x[n]} and the > Fourier series of x(t)?You observe that when you sample x(t) to get x[n], aliasing happens. Let Ts be the sample interval, N be the number of points in x[n], Xd = DFT {x[n]} and Xc = F{x(t)}. Then (possibly with some difficulties in scaling, and a very strong need to pay attention to the relative phase of x(t) and the sampling clock), taking aliasing into account gives you Xd[k] = sum_{p = -\infty}^{\infty}{Xc(2 * \pi * (p + k/N)/Ts)} Some hard thinking about the sampling process will reveal all of this to you. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
Reply by ●August 6, 20142014-08-06
Hi Tim, Thanks - first of all, do you agree with what I wrote? If not, could you point me to the specific word? Believe me, I spent hours in thinking about it. I don't see how this could be correct. So I pinpointed it by exactly reproducing word-by-word what the definiton says. On 2014-08-06 12:56, Tim Wescott wrote:> On Wed, 06 Aug 2014 12:31:01 -0700, Peter Mairhofer wrote: >> Half a year I posted here about an inconsistency between DFT and its >> continuous repeated, interpolated Fourier series which should match (up >> to a factor). I did not really get an answer to it. So I try it again - >> step by step. >> >> My main point is that this is not true except we do NOT use the >> interplated signal but rather the infinite-bandwidth, non-interpolated >> signal! >> >> So I quote the relevant sentences from >> https://ccrma.stanford.edu/~jos/mdft/Relation_DFT_Fourier_Series.html >> (or >> http://www.dsprelated.com/dspbooks/mdft/ > Relation_DFT_Fourier_Series.html) >> and use "JOS>> " as ident. >> >> JOS>> We now show that the DFT of a sampled signal x(n) (of length N) >> JOS>> ... >> >> Ok, I plotted one: x = [ 2 1 ], N=2: >> >> http://snag.gy/Q4AWO.jpg >> >> Let's call it "x[n]". >> >> JOS>> ... is proportional to the Fourier series coefficients of the >> JOS>> continuous periodic signal obtained by repeating ... >> >> Ok, converted the signal to a continuous signal and repeated: >> >> http://snag.gy/woEyc.jpg >> >> Let's call it "xs(t)". >> >> Note that: >> 1.) I zoomed into an arbitrary time window (it's periodic) >> 2.) The arrows represent weighted Diracs 3.) The (continuous) values >> between the Diracs are zero(!) >> 4.) The signal has infinite bandwidth 5.) The period T0 = 2 >> >> JOS>> ... and interpolating x. >> >> Ok, interpolated: >> >> http://snag.gy/TsfPr.jpg >> >> Let's call it "x(t)". >> >> For better clarity, the non-interpolated xs(t) is also included. >> >> Note that the interpolation corresponds to the convolution with an >> appropriately weighted sinc or equivalently, filtering with the >> appropriate box lowpass filter. >> >> *Now* the analog signal is bandlimited. >> And it is *nonzero* everywhere! >> >> And according to the quoted sentence (I read it about thousand times) >> the Fourier coefficients obtained from x(t) over one period should be >> proportional to the DFT coefficients of x[n]. >> >> Then he does a bunch of stuff, but then he writes: >> >> JOS>> We wish to find the continuous-time Fourier series of the JOS>> >> sampled periodic signal x(nT) >> >> NO!! You told something else above! We wish to find the "continuous-time >> Fourier series of the periodic, *interpolated* signal x(t)" >> >> If I continue this way - of course, the result is obvious because I >> multiply with the Dirac comb and I am back to xs(t). But that's not my >> bandlimited, analog version of x[n]. I delete the non-zero data which >> was introduced by the interpolation. >> >> I know that the relation is conceptionally correct but I insist that the >> quoted document is not written consistently (if you do not agree could >> you show why it is correct?) >> >> In any case, how can I find the relation between DFT{x[n]} and the >> Fourier series of x(t)? > > You observe that when you sample x(t) to get x[n], aliasing happens.Aehm, no? x(t) is bandlimited to [-piT,piT] so no aliasing should occur. But even if not, the question is still up: "how can I find the relation between DFT{x[n]} and the Fourier series of x(t)?" (and not xc(t)). An alternative way to think about it: You send x[n] through a DAC and filter it. What you should get is x(t) (and not xc(t)). Now do the multiplication and integration in analog - just by definition - without zeroing out anything.> Let Ts be the sample interval, N be the number of points in x[n], Xd = DFT > {x[n]} and Xc = F{x(t)}.infty x(t) = SUM x[mod(n,N)] sinc(t/Ts-n) n=-infty right? And T0 F{x(t)} = INT x(t) exp(-i*w0*k*t) dt 0 where T0 is the period and w0=2pi/T0.> Then (possibly with some difficulties in > scaling, and a very strong need to pay attention to the relative phase of > x(t) and the sampling clock), taking aliasing into account gives you > > Xd[k] = sum_{p = -\infty}^{\infty}{Xc(2 * \pi * (p + k/N)/Ts)}Yes, I think that's clear to me. But I don't see how to reconcile that with the rest. On the other hand, if I execute F{x(t)} I do not get this. I get it only if I execute F{xc(t)} where xc(t) is not bandlimited and infty xc(t) = SUM x[mod(n,N)] \delta(t/Ts-n) n=-infty ... which is what I show already in my plots.> Some hard thinking about the sampling process will reveal all of this to > you.It did not. That's why I'm desperately seeking for help. Thanks Peter
Reply by ●August 6, 20142014-08-06
On 8/6/14 4:25 PM, Peter Mairhofer wrote:> >> Some hard thinking about the sampling process will reveal all of this to >> you. > > It did not. That's why I'm desperately seeking for help. >there are two or three processes that you should think about (and it need not be all that hard, since, from your notation, you seem to get a handle on the mathematical expression). 1. sampling process. what happens to the spectrum of x(t) when you sample it with a "dirac comb" or whatever you want to call the sampling function? 2. windowing to finite length. what happens when you window either x(t) or the sampled version, x[n], with a rectangular window of length N? 3. periodic extension. what happens when you periodically extend it by repeatedly shifting x[n] by N samples and overlap and add it? deal with each step by itself. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●August 6, 20142014-08-06
On Wed, 06 Aug 2014 12:31:01 -0700, Peter Mairhofer wrote:> Hi all, > > Half a year I posted here about an inconsistency between DFT and its > continuous repeated, interpolated Fourier series which should match (up > to a factor). I did not really get an answer to it. So I try it again - > step by step. > > My main point is that this is not true except we do NOT use the > interplated signal but rather the infinite-bandwidth, non-interpolated > signal! > > So I quote the relevant sentences from > https://ccrma.stanford.edu/~jos/mdft/Relation_DFT_Fourier_Series.html > (or > http://www.dsprelated.com/dspbooks/mdft/Relation_DFT_Fourier_Series.html)> and use "JOS>> " as ident. > > JOS>> We now show that the DFT of a sampled signal x(n) (of length N) > JOS>> ... > > Ok, I plotted one: x = [ 2 1 ], N=2: > > http://snag.gy/Q4AWO.jpg > > Let's call it "x[n]". > > JOS>> ... is proportional to the Fourier series coefficients of the > JOS>> continuous periodic signal obtained by repeating ... > > Ok, converted the signal to a continuous signal and repeated: > > http://snag.gy/woEyc.jpg > > Let's call it "xs(t)". > > Note that: > 1.) I zoomed into an arbitrary time window (it's periodic) > 2.) The arrows represent weighted Diracs 3.) The (continuous) values > between the Diracs are zero(!) > 4.) The signal has infinite bandwidth 5.) The period T0 = 2 > > JOS>> ... and interpolating x. > > Ok, interpolated: > > http://snag.gy/TsfPr.jpg > > Let's call it "x(t)". > > For better clarity, the non-interpolated xs(t) is also included. > > Note that the interpolation corresponds to the convolution with an > appropriately weighted sinc or equivalently, filtering with the > appropriate box lowpass filter. > > *Now* the analog signal is bandlimited. > And it is *nonzero* everywhere! > > And according to the quoted sentence (I read it about thousand times) > the Fourier coefficients obtained from x(t) over one period should be > proportional to the DFT coefficients of x[n]. > > Then he does a bunch of stuff, but then he writes: > > JOS>> We wish to find the continuous-time Fourier series of the JOS>> > sampled periodic signal x(nT) > > NO!! You told something else above! We wish to find the "continuous-time > Fourier series of the periodic, *interpolated* signal x(t)" > > If I continue this way - of course, the result is obvious because I > multiply with the Dirac comb and I am back to xs(t). But that's not my > bandlimited, analog version of x[n]. I delete the non-zero data which > was introduced by the interpolation. > > I know that the relation is conceptionally correct but I insist that the > quoted document is not written consistently (if you do not agree could > you show why it is correct?) > > In any case, how can I find the relation between DFT{x[n]} and the > Fourier series of x(t)? > > Thanks Peter"Digital Signal Processing" by Oppenheim and Shaffer. "Signals and Systems" by Oppenheim, Wilsky with Young. They're ancient (over 30 years old), but the math hasn't changed. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
Reply by ●August 6, 20142014-08-06
Hi Robert, On 2014-08-06 13:58, robert bristow-johnson wrote:> On 8/6/14 4:25 PM, Peter Mairhofer wrote: >> >>> Some hard thinking about the sampling process will reveal all of this to >>> you. >> >> It did not. That's why I'm desperately seeking for help. >> > > there are two or three processes that you should think about (and it > need not be all that hard, since, from your notation, you seem to get a > handle on the mathematical expression).Why don't you look what I have written? I really tried to make it as crisp and clear is it can get but I feel my post is not read and hence misunterstood. Is this x[n]? http://snag.gy/Q4AWO.jpg Is this xs(t)? http://snag.gy/woEyc.jpg Is this x(t)? http://snag.gy/TsfPr.jpg (the green one) JOS wrote: "We now show that the DFT of a sampled signal x(n) (of length N) is proportional to the Fourier series coefficients of the continuous periodic signal obtained by repeating and interpolating x." Is he talking about x(t)? Then there is an inconsistency ... somewhere. It might be that there is a different way to look at it but then I do not understand WHY you obtain different results when just using plain definitions from time domain signals. However, overall I think that we are talking about different things. I do not want to sample. I want to obtain the Fourier series of an analog signal. Without sampling. Nothing else was written. x(t) = "DAC(x[n]) + LOWPASS" Integral of x(t) and NOT x(nT).> 1. sampling process. what happens to the spectrum of x(t) when you > sample it with a "dirac comb" or whatever you want to call the sampling > function?As I said, I do not want to sample. But if I did, I get exactly x[n] because infty x(t) = SUM x[n] sinc(t/Ts-n) n=-infty Yes, if I sample and then use x(nTs) the result is obvious. I wrote that.> 2. windowing to finite length. what happens when you window either > x(t) or the sampled version, x[n], with a rectangular window of length N?What happens to what? The signal has length N and/or spectrum gets leakage (sinc convolution)> 3. periodic extension. what happens when you periodically extend it > by repeatedly shifting x[n] by N samples and overlap and add it?Zero padding in frequency domain. But I do not see (yet) how this relates to the question or how to reconcile it. In the end it's a trivial yes/no question. Is the green signal here http://snag.gy/TsfPr.jpg (which is x(t)) what you/Julius Smith is referring to to integrate per definition of the the Fourier series integral? If no, we do not yet talk about the same thing. If yes, the text is at least inconsistent. Thanks Peter
Reply by ●August 6, 20142014-08-06
On 2014-08-06 14:22, Tim Wescott wrote:> On Wed, 06 Aug 2014 12:31:01 -0700, Peter Mairhofer wrote: > [...] > "Digital Signal Processing" by Oppenheim and Shaffer.^^^^^^^ Schafer. Have it, read it - already 9 years ago. I do not expect anyone answering here and I really appreciate when people spend time in helping others (you provided already very helpful answers in other topics). On the other hand, when you do not read my posting and/or not refer to its content then you do not only not help me but also waste your time - which is also valuable I guess. Peter PS: OT, hence f'up 2 poster
Reply by ●August 6, 20142014-08-06
Hi Robert, On 2014-08-06 13:58, robert bristow-johnson wrote:> On 8/6/14 4:25 PM, Peter Mairhofer wrote: >> >>> Some hard thinking about the sampling process will reveal all of this to >>> you. >> >> It did not. That's why I'm desperately seeking for help. >> > > there are two or three processes that you should think about (and it > need not be all that hard, since, from your notation, you seem to get a > handle on the mathematical expression). > > 1. sampling process. what happens to the spectrum of x(t) when you > sample it with a "dirac comb" or whatever you want to call the sampling > function? > > 2. windowing to finite length. what happens when you window either > x(t) or the sampled version, x[n], with a rectangular window of length N? > > 3. periodic extension. what happens when you periodically extend it > by repeatedly shifting x[n] by N samples and overlap and add it? > > deal with each step by itself.I think I know what you want to tell me here - please let me know if this is the case: Start off with a continous signal x(t) which is periodic with T0 and bandlimited to 1/2Ts. Its Fourier series is T0 Xk = integral x(t) exp(-j2pi/T0 k t) dt 0 If I would sample it at fs=1/Ts, that corresponds to a multiplication with the Dirac comb which is xs(t) from my posting. Plugging this in the above definition results in T0 Xk = integral xc(t) exp(-j2pi/T0 k t) dt 0 T0 N-1 Xk = integral SUM x(nTs) \delta(t-nTs) exp(-j2pi/T0 k t) dt 0 n=0 N-1 T0 = SUM x(nTs) integral \delta(t-nTs) exp(-j2pi/T0 k t) dt n=0 0 \-----------------------------------------/ this is only nonzero when t=nTs so we can set t=nTs and x(nTs) = x[n] N-1 = SUM x[n] exp(-j2pi/N kn) = DFT{x[n]} n=0 This is clear. But the confusion comes when I go the opposite direction. I start with x[n] that has a certain DFT. From this I create first the periodic extension xp[n]=x[mod(n,N)], n=-infty,...infty, from this the analog signal infty xs(t) = SUM xp[n] \delta(t-nTs) n=-infty and then the actual signal x(t) by low-passing it. infty x(t) = SUM xp[n] sinc(t/Ts-n) n=-infty Why, when I now execute the definition of the continuous Fourier series on x(t) do I not get back to DFT{x[n]} again? (as mentioned multiple times, I do if I would execute it on xs(t) ) Regards, Peter
Reply by ●August 6, 20142014-08-06
On Wed, 06 Aug 2014 15:28:38 -0700, Peter Mairhofer wrote:> Hi Robert, > > On 2014-08-06 13:58, robert bristow-johnson wrote: >> On 8/6/14 4:25 PM, Peter Mairhofer wrote: >>> >>>> Some hard thinking about the sampling process will reveal all of this >>>> to you. >>> >>> It did not. That's why I'm desperately seeking for help. >>> >>> >> there are two or three processes that you should think about (and it >> need not be all that hard, since, from your notation, you seem to get a >> handle on the mathematical expression). >> >> 1. sampling process. what happens to the spectrum of x(t) when >> you >> sample it with a "dirac comb" or whatever you want to call the sampling >> function? >> >> 2. windowing to finite length. what happens when you window either >> x(t) or the sampled version, x[n], with a rectangular window of length >> N? >> >> 3. periodic extension. what happens when you periodically extend >> it >> by repeatedly shifting x[n] by N samples and overlap and add it? >> >> deal with each step by itself. > > I think I know what you want to tell me here - please let me know if > this is the case: > > Start off with a continous signal x(t) which is periodic with T0 and > bandlimited to 1/2Ts. > > Its Fourier series is > > T0 > Xk = integral x(t) exp(-j2pi/T0 k t) dt > 0 > > If I would sample it at fs=1/Ts, that corresponds to a multiplication > with the Dirac comb which is xs(t) from my posting. Plugging this in the > above definition results in > > T0 > Xk = integral xc(t) exp(-j2pi/T0 k t) dt > 0 > > > T0 N-1 > Xk = integral SUM x(nTs) \delta(t-nTs) exp(-j2pi/T0 k t) dt > 0 n=0 > > N-1 T0 > = SUM x(nTs) integral \delta(t-nTs) exp(-j2pi/T0 k t) dt > n=0 0 > \-----------------------------------------/ > this is only nonzero when t=nTs > so we can set t=nTs > > and x(nTs) = x[n] > > N-1 > = SUM x[n] exp(-j2pi/N kn) = DFT{x[n]} > n=0 > > This is clear. > > But the confusion comes when I go the opposite direction. I start with > x[n] that has a certain DFT. From this I create first the periodic > extension xp[n]=x[mod(n,N)], n=-infty,...infty, from this the analog > signal > > infty > xs(t) = SUM xp[n] \delta(t-nTs) > n=-infty > > and then the actual signal x(t) by low-passing it. > > infty > x(t) = SUM xp[n] sinc(t/Ts-n) > n=-infty > > Why, when I now execute the definition of the continuous Fourier series > on x(t) do I not get back to DFT{x[n]} again? > > (as mentioned multiple times, I do if I would execute it on xs(t) )You should, indeed, get a zero-padded and possibly scaled version of DFT{x[n]}. When I get into this sort of bind I go all the way back to the basics. Try it with one cycle of a sine wave. Try it with a known 1st- and 2nd- harmonic added together, etc. The theory, as stated, should work, so there is something in your execution that is creating difficulty. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
Reply by ●August 7, 20142014-08-07
On 2014-08-06 16:34, Tim Wescott wrote:> On Wed, 06 Aug 2014 15:28:38 -0700, Peter Mairhofer wrote: >> [...] >> >> Start off with a continous signal x(t) which is periodic with T0 and >> bandlimited to 1/2Ts. >> >> Its Fourier series is >> >> T0 >> Xk = integral x(t) exp(-j2pi/T0 k t) dt >> 0 >> >> If I would sample it at fs=1/Ts, that corresponds to a multiplication >> with the Dirac comb which is xs(t) from my posting. Plugging this in the >> above definition results in >> >> T0 >> Xk = integral xc(t) exp(-j2pi/T0 k t) dt >> 0 >> >> >> T0 N-1 >> Xk = integral SUM x(nTs) \delta(t-nTs) exp(-j2pi/T0 k t) dt >> 0 n=0 >> >> N-1 T0 >> = SUM x(nTs) integral \delta(t-nTs) exp(-j2pi/T0 k t) dt >> n=0 0 >> \-----------------------------------------/ >> this is only nonzero when t=nTs >> so we can set t=nTs >> >> and x(nTs) = x[n] >> >> N-1 >> = SUM x[n] exp(-j2pi/N kn) = DFT{x[n]} >> n=0 >> >> This is clear. >> >> But the confusion comes when I go the opposite direction. I start with >> x[n] that has a certain DFT. From this I create first the periodic >> extension xp[n]=x[mod(n,N)], n=-infty,...infty, from this the analog >> signal >> >> infty >> xs(t) = SUM xp[n] \delta(t-nTs) >> n=-infty >> >> and then the actual signal x(t) by low-passing it. >> >> infty >> x(t) = SUM xp[n] sinc(t/Ts-n) >> n=-infty >> >> Why, when I now execute the definition of the continuous Fourier series >> on x(t) do I not get back to DFT{x[n]} again? >> >> (as mentioned multiple times, I do if I would execute it on xs(t) ) > > You should, indeed, get a zero-padded and possibly scaled version of > DFT{x[n]}. > > When I get into this sort of bind I go all the way back to the basics. > Try it with one cycle of a sine wave. Try it with a known 1st- and 2nd- > harmonic added together, etc. The theory, as stated, should work, so > there is something in your execution that is creating difficulty.Thanks Tim; I did that. Finally, after hours I got it. The main thing was the ordering of the coefficients - I forgot to realize that k goes from -floor(K/2)+1 to floor(K/2) rather than 0,...,K-1. The minor thing is that I started with "N=1" and always got correct results from the beginning on. This is obvious because it's just DC (the cos evaluates to 1, the sine to 0). Then I continued with N=2 and got different results than expected. The weird thing now is: It works for all N except for N=2. What could possibly the cause for that? I show a small piece of MATLAB code here: N = 3; xn = ifft([ 7 2 3 ]'); DFTx = fft(xn); Periods = 100; xp = repmat(xn,Periods,1); % analog, sampled signal Ts = 1; fs = 1/Ts; SimFactor = 1000; xs = upsample(xp,SimFactor); xt = resample(xp,SimFactor,1); % "simulation" sampling rate, would go -> 0 in real setup Tss = Ts/SimFactor; T0 = Ts*N; t = (0:Tss:Ts*N*Periods-Tss)'; % integration interval, pick one period from the middle to % minimize edge effects ta = floor(Periods/2)*T0; tb = ta + T0; idx = find(t>=ta & t<tb); FSxt = []; for k=[ -1 0 1 ] % evaluate analog definition of Fourier series xmod = xt .* exp(-1i*2*pi/T0*k*t); dt = 1/SimFactor; FSxt = [FSxt ; 1/T0 * sum(xmod(idx)) * dt ]; end [ DFTx , N*FSxt ] This is the result: [ DFTx , N*FSxt ] = 7.0000 6.9953 - 0.0000i 2.0000 2.0022 - 0.0000i 3.0000 3.0033 + 0.0000i which differs by ~0.1% but I still assume a match and attribute the discrepancy to simulation artefacts. If I change N to 1 or any number >2 it works too. However, when I set N=2 and replace xn = ifft([ 7 2 3 ]'); -> xn = ifft([ 7 2 ]'); for k=[ 0 1 -1 ] -> for k=[ 0 1 ] the second bin is always HALF of what it should be, e.g.: [ DFTx - N*FSxt ] = 7.0000 6.9953 2.0000 0.9996 + 0.0000i It must be something trivial but since it's only N=2 I don't see what! Furthermore it still looks a little bit "magic" that it works. I mentioned a couple of times, I immideately see the correct result for "xs(t)" but for not for x(t). I tried to evaluate the integral over x(t) as defined at the top of this posting but sure enough the integral does not even converge (it's essentially the integral over the sinc). So why does it indeed work (except for N=2) and why does integrating over x(t) fail? (And my original question is up too, why Julius Smith talks about x(t) but entirely operates on xs(t) which is still contradicting to me) Thanks, Peter
