So, this was an ongoing discussion in another thread, and it's been going
on long enough that it probably deserves its own thread. I'm going to
switch the notation to use 'n' in the time domain and 'k' in the frequency
domain.
To summarize the problem:
Given a signal x[n] such that x[n+N] = x[n] for all n. x[n] may be
complex.
Provide a continuous expression x(t) such that x(t) = x[n] for t = n,
that is bandwidth limited by the Nyquist frequency.
The method I suggested uses a DFT. The idea is to use the IDFT to obtain
an expression for x[n] that uses sines and cosines.
So the DFT shall be defined as,
X[k] = sum( x[n]*e^(-j*2pi*n*k/N) ) for n = 0 to N-1 (1)
The IDFT shall be defined as,
x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = 0 to N-1 (2)
Now,
e^(-j*2pi*n*(k+N)/N) = e^(-j*2pi*n*k/N) * e^(-j*2pi*n)
Because n is always an integer value,
e^(-j*2pi*n*k/N) = e^(-j*2pi*n*(k+N)/N) (3)
So by extension,
X[k] = X[k+N] (4)
This allows us to modify our definition for the IDFT using (3) and (4),
x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = c to c+N-1 (5)
Where c is any integer. For the purposes of this exercise, it is ideal to
pick c such that all values of k fall within the Nyquist bandwidth. Lets
start with the case where N is odd, then c = -(N-1)/2. Rewriting (5),
x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = -(N-1)/2 to (N-1)/2
(6)
We now have an expression for x[n] using sines and cosines that can be made
continuous,
x(t) = (1/N) * sum( X[k]*( cos(2pi*t*k/N) + j*sin(2pi*t*k/N) ) ) for k
= -(N-1)/2 to (N-1)/2 (7)
Lets examine the case where x[n] is real. In this case,
X[-k] = X*[k] (8)
Where X*[k] is the complex conjugate of X[k]. To simplify this expression,
lets look at sum of a given k and -k.
X[k] * ( cos(2pi*t*k/N) + j*sin(2pi*t*k/N) ) + X[-k] * (
cos(-2pi*t*k/N) + j*sin(-2pi*t*k/N) )
= 2 * Re(X[k]) * cos(2pi*t*k/N) - 2 * Im(X[k]) * sin(2pi*t*k/N) (9)
Now we can factor out X[0] from (7) and be left with pairs of positive and
negative k's that can be factored using (9),
x(t) = (1/N) * ( X[0] + 2*sum( Re(X[k])*cos(2pi*t*k/N) -
Im(X[k])*sin(2pi*t*k/N) ) ) for k = 1 to (N-1)/2 (10)
Equation (10) is the solution for the case where x[n] is real and N samples
long where N is odd. The derivation for the case where N is even is
similar. Lets start with (5) with c = -N/2.
x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = -N/2 to (N/2)-1
(11)
This need to made to have k pairs at all positive and negative frequencies.
From (4),
X[-N/2] = (1/2) * ( X[-N/2] + X[N/2] ) (12)
Plugging (12) into (11) and making it continuous,
x(t) = (1/N) * ( X[N/2]*cos(pi*t) + sum( X[k]*(cos(2pi*t*k/N) +
j*sin(2pi*t*k/N)) ) ) for k = -((N/2)-1) to (N/2)-1 (13)
This can be applied for a real x[n] in a similar fashion as (10)
x(t) = (1/N) * ( X[0] + X[N/2]*cos(pi*t) + 2*sum(
Re(X[k])*cos(2pi*t*k/N) - Im(X[k])*sin(2pi*t*k/N) ) ) for k = 1 to (N/2)-1
(14)
That about does it... Equations (7) and (13) can be used for the case
where x[n] is complex, and reduce to equations (10) and (14) where x[n] is
real.
I welcome any feedback on the usefulness, cleanliness and/or clarity of
this post.
Cheers,
-Dan
Interpolation of a Discrete Periodic Signal
Started by ●May 23, 2013
Reply by ●May 23, 20132013-05-23
On May 23, 12:05�am, "dszabo" <62466@dsprelated> wrote:> So, this was an ongoing discussion in another thread, and it's been going > on long enough that it probably deserves its own thread. �I'm going to > switch the notation to use 'n' in the time domain and 'k' in the frequency > domain. > > To summarize the problem: > > � � Given a signal x[n] such that x[n+N] = x[n] for all n. �x[n] may be > complex. > > � � Provide a continuous expression x(t) such that x(t) = x[n] for t = n, > that is bandwidth limited by the Nyquist frequency. > > The method I suggested uses a DFT. �The idea is to use the IDFT to obtain > an expression for x[n] that uses sines and cosines. > > So the DFT shall be defined as, > > � � X[k] = sum( x[n]*e^(-j*2pi*n*k/N) ) for n = 0 to N-1 � �(1) > > The IDFT shall be defined as, > > � � x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = 0 to N-1 � �(2) > > Now, > > � � e^(-j*2pi*n*(k+N)/N) = e^(-j*2pi*n*k/N) * e^(-j*2pi*n) > > Because n is always an integer value, > > � � e^(-j*2pi*n*k/N) = e^(-j*2pi*n*(k+N)/N) � �(3) > > So by extension, > > � � X[k] = X[k+N] � �(4) > > This allows us to modify our definition for the IDFT using (3) and (4), > > � � x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = c to c+N-1 � �(5) > > Where c is any integer. �For the purposes of this exercise, it is ideal to > pick c such that all values of k fall within the Nyquist bandwidth. �Lets > start with the case where N is odd, then c = -(N-1)/2. �Rewriting (5), > > � � x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = -(N-1)/2 to (N-1)/2 > � (6) > > We now have an expression for x[n] using sines and cosines that can be made > continuous, > > � � x(t) = (1/N) * sum( X[k]*( cos(2pi*t*k/N) + j*sin(2pi*t*k/N) ) ) for k > = -(N-1)/2 to (N-1)/2 � �(7) > > Lets examine the case where x[n] is real. �In this case, > > � � X[-k] = X*[k] � �(8) > > Where X*[k] is the complex conjugate of X[k]. �To simplify this expression, > lets look at sum of a given k and -k. > > � � X[k] * ( cos(2pi*t*k/N) + j*sin(2pi*t*k/N) ) + X[-k] * ( > cos(-2pi*t*k/N) + j*sin(-2pi*t*k/N) ) > > � � = 2 * Re(X[k]) * cos(2pi*t*k/N) - 2 * Im(X[k]) * sin(2pi*t*k/N) � �(9) > > Now we can factor out X[0] from (7) and be left with pairs of positive and > negative k's that can be factored using (9), > > � � x(t) = (1/N) * ( X[0] + 2*sum( Re(X[k])*cos(2pi*t*k/N) - > Im(X[k])*sin(2pi*t*k/N) ) ) for k = 1 to (N-1)/2 � �(10) > > Equation (10) is the solution for the case where x[n] is real and N samples > long where N is odd. �The derivation for the case where N is even is > similar. �Lets start with (5) with c = -N/2. > > � � x[n] = (1/N) * sum( X[k]*e^(j*2pi*n*k/N) ) for k = -N/2 to (N/2)-1 > (11) > > This need to made to have k pairs at all positive and negative frequencies. > �From (4), > > � � X[-N/2] = (1/2) * ( X[-N/2] + X[N/2] ) � �(12) > > Plugging (12) into (11) and making it continuous, > > � � x(t) = (1/N) * ( X[N/2]*cos(pi*t) + sum( X[k]*(cos(2pi*t*k/N) + > j*sin(2pi*t*k/N)) ) ) �for k = -((N/2)-1) �to (N/2)-1 � �(13) > > This can be applied for a real x[n] in a similar fashion as (10) > > � � x(t) = (1/N) * ( X[0] + X[N/2]*cos(pi*t) + 2*sum( > Re(X[k])*cos(2pi*t*k/N) - Im(X[k])*sin(2pi*t*k/N) ) ) for k = 1 to (N/2)-1 > � (14) > > That about does it... �Equations (7) and (13) can be used for the case > where x[n] is complex, and reduce to equations (10) and (14) where x[n] is > real. > > I welcome any feedback on the usefulness, cleanliness and/or clarity of > this post. >Dan, i think it's all accurate. one thing that it looks like you're saying (with Eq. 4) is that the DFT is inherently about periodic discrete sequences. you hadn't said that x[n+N]=x[n] like X[k +N]=X[k], but it's true and pushing this fact has gotten me into tiffs here on this newsgroup. i, essentially, refuse to differentiate between the concepts of the DFS and the DFT. now, one thing you can do to "clean" it a little is to make use of the finite geometric series: K-1 SUM{ p^k } = (p^K - 1)/(p-1) for any p k=0 and apply that to the simpler case where x[n] = {1, 0, 0, 0, ...0} for 0 <= n < N (as you have previously noted, all of the X[k]=1). and that simpler case can be extended to a general periodic x[n] using linearity and time-invariancy. if you do that, you can get something that looks exactly like the Dirichlet kernel for N odd: N-1 x(t) = SUM{ x[n] sin(pi*(t-n)) / ( N*sin((pi/N)*(t-n)) ) } n=0 there is a problem with the N even case because there is one component of X[k] that is precisely on the Nyquist frequency, namely X[N/2]. now, if you do all of that geometric summation crap, you have to decide what to do with the X[N/2] component. leave it in the positive frequencies or put it in the negative frequencies? well, if all x[n] are real, we know that X[N-k] = X*[k] which means X[0] must be real and X[N/2], if it exists, is real. but the bandlimited x(t) *could* have an imaginary component to that (which is opposite sign for -N/2 than for N/2) and when x(t) is hypothetically sampled to be x[n], the spectrum is overlap-added and that imaginary component at -N/2 will be cancelled by the imaginary component at +N/ 2. but the real components *must* add up. so then, for the continuous x(t), half of X[N/2] must go down to X[-N/2]. if you leave the imaginary part as zero you get as a result, for N even: N-1 x(t) = SUM{ x[n] sin(pi*(t-n)) / ( N*tan((pi/N)*(t-n)) ) } n=0 notice the tan() in the denominator instead of sin(). so it's a little different from that Dirichlet thingie. if anyone insists on an arbitrary imaginary component to the Nyquist component, they will get: N-1 x(t) = SUM{ x[n] sin(pi*(t-n)) / ( N*tan((pi/N)*(t-n)) ) } + A*sin(pi*t) n=0 A can be any number and it will still add to x[n] when t=n. but if x(t) is linearly dependent on x[n], then A *must* be proportional to the sum of all of the x[n] so that when x[n] goes to zero, so also will A (or when all of the x[n] is doubled, so also does A). but i think that it's best for A to be zero. in either the N even or N odd cases, it's like we have little periodic sinc() functions in the interpolation, as would be expected. for the moment, i am stuck with Google Groups and i cannot predict when they insert the 3D in with the = sign or not, so if this is unreadable with your newsreader, try viewing with GG. r b-j
Reply by ●May 23, 20132013-05-23
Ohhhh yeah, totally. I hadn't though of trying to reduce it to using the actual signal terms to build it. It reminds me of converting a moving average to a CIC filter. I'm gonna try and look at that this weekend if I have time!
Reply by ●May 24, 20132013-05-24
On Thursday, May 23, 2013 3:29:24 PM UTC-4, dszabo wrote:> I hadn't though of trying to reduce it to using the > actual signal terms to build it.it was the original problem statement. find the function: x(t) = function({x[0],x[1],x[2],..x[N-1]}, t)> It reminds me of converting a moving average to a CIC filter.i always thought a moving average *was* a CIC filter. do you mean converting the "naive" transversal structure (with all these coefficients that are 1/L) into the recursive structure (that has less computation)? r b-j
Reply by ●June 3, 20132013-06-03
Don't know if it originated as a homework problem, but it would make a good one. Anwyay, the word is out. Simply write out the equations of the IFFT with continuous time replacing sample indices, and you can interpolate the signal anywhere, with arbitrary precision. Showing that it is "the" bandlimited signal - Nyquist, there can be only one - that crosses through all the sampled points Or, "almost". It works fine when n is odd. But, the final bin in an even-sized FFT doesn't play along. Figuring this out is a bonus question, and a nasty one.
Reply by ●June 3, 20132013-06-03
On Jun 3, 3:30�pm, "mnentwig" <24789@dsprelated> wrote:> Don't know if it originated as a homework problem, but it would make a good > one. Anwyay, the word is out. > > Simply write out the equations of the IFFT with continuous time replacing > sample indices, and you can interpolate the signal anywhere, with arbitrary > precision. > Showing that it is "the" bandlimited signal - Nyquist, there can be only > one - that crosses through all the sampled points > > Or, "almost". It works fine when n is odd. But, the final bin in an > even-sized FFT doesn't play along. Figuring this out is a bonus question, > and a nasty one.aw, split the Nyquist bin, X[N/2], in half. send half of it to the negative frequency and and keep the remaining half at the positive frequency. the real part (which results in A*cos(pi*n)) is well determined. A is known. the imaginary part (which results in B*sin(pi*n)) is completely indeterminate. B can be any number and the interpolated function will still go through all of the sampled points. but i would set B to zero. there are a few other "homework problems" i can suggest. even MSEE problems. r b-j
Reply by ●June 4, 20132013-06-04
>aw, split the Nyquist bin, X[N/2], in half. send half of it to the >negative frequency and and keep the remaining half at the positive >frequency.That would be a funny OFDM spectrum if I start with a vector full of subcarriers :-) Problem is, we don't know whether it's positive or negative frequency or both, this needs (in my interpretation) extra information.
Reply by ●June 4, 20132013-06-04
On Jun 4, 1:46�am, "mnentwig" <24789@dsprelated> wrote:> >aw, split the Nyquist bin, X[N/2], in half. �send half of it to the > >negative frequency and and keep the remaining half at the positive > >frequency. > > That would be a funny OFDM spectrum if I start with a vector full of > subcarriers :-)i was originally restricting this periodic interpolation to real x[n] and real x(t).> Problem is, we don't know whether it's positive or negative frequency or > both, this needs (in my interpretation) extra information.if x(t) is real, we know. it's both and the positive and negative frequency components have hermitian symmetry (are complex conjugates of each other). think of it in reverse for a minute: what must X(w) be if x(t) is real and is bandlimited to |w| <= pi? (let's assume that the sampling period is 1.) when x(t) is sampled at integer times (t=n), then what happens to the spectrum? this doesn't settle what the amplitude of the sin(pi*n) Nyquist component could be (it could be anything and it would not change the value of x(t) when t=n), but i am still happy to make that amplitude zero by fiat. however we know what the cos(pi*n) component must be. we have enough information for that (if x(t) is real). r b-j
Reply by ●June 4, 20132013-06-04
Yes, but it's different, depending whether you get the original FFT vector from transforming (and implicitly sampling) a real signal or a complex one. For a complex one, I can't distinguish after sampling, whether the highest frequency bin is spinning forwards or backwards. So I have to know, whether the signal was was bandlimited to -fs/2 <= f < fs/2 or -fs/2 < f <= fs/2 or -fs/2 <= f <= fs/2 knowing it's real-valued. This seems like splitting hairs, but it makes the difference, whether my simulator noise floor is, say, -90 dB or -300 dB.
Reply by ●June 4, 20132013-06-04
On Jun 4, 11:01�am, "mnentwig" <24789@dsprelated> wrote:> Yes, but it's different, depending whether you get the original FFT vector > from transforming (and implicitly sampling) a real signal or a complex one. > > For a complex one, I can't distinguish after sampling, whether the highest > frequency bin is spinning forwards or backwards.and it doesn't make any difference *at*the*sampling*instances* (but in between it does).> So I have to know, whether > the signal was was bandlimited to > -fs/2 <= f < fs/2 > or > -fs/2 < f <= fs/2 > or > -fs/2 <= f <= fs/2 knowing it's real-valued.i am certain that this works, because even if x[n] is complex, the interpolation kernel should be linear and time-invariant. that means, for a periodic sequence of period N, you can start with this prototype function and build offa it: { 1 for n = m*N where m = any integer p[n] = { { 0 otherwise the DFT (or DFS, it's the same thing even if Eric and Dale disagree) of that periodic sequence is P[k] = 1 for all k. if you want to do bandlimited reconstruction of x[n] (we'll call it x(t) to show it as a continuous-time function -> x(t) = x[n] when t = n = integer) then, whether x[n] is real or complex, you can build up a general x(t) from scaled and shifted versions p(t) +inf p(t) = SUM{ sinc(t - m*N) } m=-inf if you expand the sin(pi*(n-m*N)) in the numerator and take advantage of some trig identities, you get (for either even or odd N): +inf p(t) = sin(pi*t) * SUM{ ((-1)^(m*N))/(pi*(t - m*N)) } m=-inf which has that sin(pi*n) factor in common with the top equations and different expressions for N odd ((-1)^(m*N) toggles) and N even ((-1)^(m*N) is always 1). now if you do the DFT thing with P[k] and set the FT of p(t) to P[k] at the appropriate frequencies then for N odd, there is no Nyquist component to screw things up and you get: p(t) = sin(pi*t)/(N*sin(pi*t/N)) N odd but for even N you get p(t) = sin(pi*t)/(N*tan(pi*t/N)) N even now this is done for that prototype p[n] to p(t) which *is* real. and i believe, in either even or odd case, that p(-t) = p(t) but using LTI, you should be able to make this to work for any periodic x[n] where: x[n+N] = x[n] for all integer n. then N-1 x[i] = SUM{ x[n] * p[i-n] (integer i) n=0 and then N-1 x(t) = SUM{ x[n] * p(t-n) } n=0 or, when you substitute p(t): N-1 x(t) = SUM{ x[n] * sin(pi*(t-n)) / ( N*sin((pi/N)*(t- n)) ) } N odd n=0 N-1 x(t) = SUM{ x[n] * sin(pi*(t-n)) / ( N*tan((pi/N)*(t- n)) ) } N even n=0 BTW, i found from equating the two result pairs we get (in a very indirect manner) the following mathematical identities: +inf 1/tan(x) = SUM{ 1/(x - m*pi) } m=-inf +inf 1/sin(x) = SUM{ ((-1)^m)/(x - m*pi) } m=-inf ain't that weird? r b-j
