4
$\begingroup$

I'm teaching myself about Fourier Series and the DFT and trying to draw a stylised $\pi$ symbol by fourier epicycles as detailed by Mathologer on youtube (from 18:39 onwards), and the excellent explanations with extraordinary animations by 3Blue1Brown on youtube.

The goal is to generate something like this: enter image description here

using complex fourier series :

$$z(t)=\sum\limits_{k=-\infty}^{\infty}{c_k \, e^{ikt}}$$

with complex coefficients:

$$c_k=\frac{1}{2\pi}\int_\limits{-\pi}^{\pi}z(t) \, e^{-ikt} \, \mathrm{d}t$$

I have been able to generate an 'embryonic' $\pi$ shape for $c_k=-2 < k < 2$ and get same result as Mathologer (@19:19) but only because he provdes the five $c_k$ values (@20:12). Here's my output: Embryonic \pi shape by fourier

So back to the objective: I've created my own 120 point coordinate set for the $\pi$ symbol:

\pi coordinates

My question is how to find all the coefficients? I think the input coordinates need to be equally spaced samples suitable for input to the DFT, but despite much searching I'm still not sure what the process is from here?

PROGRESS UPDATE #3:

I've had a field day, made heaps of progress in MATLAB on the various algorithms. To distinguish output from the input $z$, I'm using $z_n$ for the $N=120$ complex sample points $z(1),z(2), ... z(N)$, and $z_t$ for the $D=180$ complex results $z_t(1),z_t(2), ... z_t(D)$ after the inverse DFT. Here's my plot for $z_t$ plus an overlay for the random point $z_t(93)$ showing the component summation arms and associated circles/epicycles (Notice the 180 points are closer together than the original 120 plotted above): pi symbol with zt(93) example

The following shows $z_t$ for $D=180$ overlaid with $z_n$ to amplify the inaccuracies, and zoomed in: pi DFT detail Still have some way to go; I really want to document the solution mathematically and experiment with ways to improve the accuracy of the resulting symbol. But I'm sensing I've crossed the mountain top, now it's just a case of toboganing all the way down! (famous last words :)

TIA for any further guidance

PS: here's a link my coordinates of sample points (since uploaded by @Olli as an Answer below, thank you sir). Each row has one $(x,y)$ pair, 120 rows:

link to ZIP file in my public dropbox folder

Here is the MATLAB program that r b-j kludged together to draw it (since updated by Chris) EVEN case first:

clearvars; close all
data = csvread("pi.csv"); % 121 rows with last repeating first
N = length(data) - 1;   % Chris added minus 1

inx = data(1:N,1);       % Chris was (:,1)
iny = data(1:N,2);       % ditto

Xk = fft(inx)/N;
Yk = fft(iny)/N;

X1 = Xk(1 : 1 + (N/2-1)     ); 
X4 = Xk(    1 + (N/2+1) : N );

% The main correction was here for X and Y: 
% the Nyquist freq must be allocated to one bin not two (previously)
Xnyq = Xk(1 + N/2);
X = [X1; Xnyq; X4];

Y1 = Yk(1 : 1 + (N/2-1)    );
Y4 = Yk(    1 + (N/2+1) : N);

Ynyq = Yk(1 + N/2); 
Y = [Y1; Ynyq; Y4];

x = N*ifft(X);
y = N*ifft(Y);

load('pi_zt_coords')
xt = real(ztt);
yt = imag(ztt);

plot(inx, iny,'o-','markersize',8)
hold on; grid on
plot(xt,yt,'k.-','markersize',8)
plot(x,  y,'mx')

xlim([100,250])
ylim([100,250])

legend('(x_{in} y_{in})','(x_t,y_t)','(x,y)','location','SouthEast')

title (['Even N =',num2str(N)]);

here is the result:

parametric even case

here is the same, but with one point removed so that NN is odd. note that there is no Nyquist value to split into two (since updated by Chris) ODD Case:

clearvars; close all

data = csvread("pi.csv");   % 121 rows with last repeating first
%data= csvread("pi_bandlimited.csv"); % from Olli's script - works too

data = vertcat(data(1:111,:), data(113:end,:)); % Delete row 112 to make N odd = 119

N = length(data) - 1;   % Chris added minus 1

inx = data(1:N,1);      % Chris (1:N,1) was (:,1)
iny = data(1:N,2);      % ditto

Xk = fft(inx)/N;
Yk = fft(iny)/N;

X1 = Xk(1 : 1 + (N-1)/2); 
X2 = Xk(1 + (N+1)/2 : N  );
X = [X1; X2];

Y1 = Yk(1 : 1 + (N-1)/2); 
Y2 = Yk(1 + (N+1)/2 : N);
Y = [Y1; Y2];

x = N*ifft(X);
y = N*ifft(Y);

load('pi_zt_coords')
xt = real(ztt);
yt = imag(ztt);

plot(inx, iny,'o-','markersize',8)
hold on; grid on
plot(xt,yt,'k.-','markersize',8)
plot(x,  y,'mx')

xlim([100,250])
ylim([100,250])

legend('(x_{in} y_{in})','(x_t,y_t)','(x,y)','location','SouthEast')

title (['Odd N = ',num2str(N)]);

and here's the result for ODD case: parametric odd case

And here's a link to the .mat file of the 180 $z_t$ coordinates: https://www.dropbox.com/s/gifbbvyfl0unv3f/pi_zt_coords.zip?dl=0

$\endgroup$
  • $\begingroup$ well, you cannot draw this as $r = f(\theta)$ because it's not a function. but you can draw this parametrically as $\Big( x(t), y(t) \Big)$ and both $x(t)$ and $y(t)$ are periodic functions of $t$. $t$ could be normalized to length of arc, that is an equal change of $t$ would result in an equal change in arc length, but i wouldn't necessarily do that. you have some sharp corners that should be explicitly anchored. however, if the arc between the sharp corners is sampled with points as equally-spaced as possible, my thinking is that the parametric equations for $x(t)$ and $y(t)$ will be as $\endgroup$ – robert bristow-johnson Jun 23 at 1:54
  • $\begingroup$ ... well behaved as possible. if you really wanna do this with polar coordinates, you can convert $x(t)$ and $y(t)$ into $r(t)$ and $\theta(t)$, but i don't know what would be advantageous about that. furthermore, if you are defining an outline font, you want to define a baseline (where $y=0$) that characters that project below (like "g", "j", "p", "q" "y") will have $y(t)<0$ for some values of $t$. also you want to define an origin point at the left so that $x(t)\ge 0$ for all $t$, unless there is character "kerning". i dunno how to deal with kerned characters. $\endgroup$ – robert bristow-johnson Jun 23 at 2:01
  • $\begingroup$ again, Chris, the points don't have to be equally-spaced (i presume you mean along the arc length) for parametric equations $\Big( x(t), y(t) \Big)$ (or for $\Big( r(t), \theta(t) \Big)$ if you prefer polar coordinates) because of the third variable $t$, which is your only "independent variable". If the points are equally-spaced, it just means that your pen is moving at a constant speed as it draws the curve. If not equally spaced, it means that the pen speeds up (where the spacing is larger) and slow down (where the spacing is less). it does not change the curve. $\endgroup$ – robert bristow-johnson Jun 23 at 4:30
  • $\begingroup$ i cannot change my last comment. i need to modify my last sentence to: "it does not change the curve where it hits the points." in-between the points the curve will be different, given different spacings of the points. $\endgroup$ – robert bristow-johnson Jun 23 at 5:10
  • 1
    $\begingroup$ @robertbristow-johnson, My followup is in response to your comments to my answer and a more thorough answer for the OP. $\endgroup$ – Cedron Dawg Jun 24 at 21:52
3
$\begingroup$

I'm not understanding the comments.

Of course you can do this. It is simply a matter of understanding what a DFT means, how to calculate DFT bin values, and how to interpret those bin values as continuous fourier series coefficients.

First off, the plane you are looking at is the complex plane. Your points are a set of $N$ discrete samples. Each sample is a complex point. Therefore what you have is the representation of one cycle of a repeating complex signal. The spacing in the diagram is not that important.

  • Any sequence of $N$ points can be represented by N coefficients exactly at the sample points

The question is: "Does your shape allow you to disregard coefficients so you have a much smaller number of them?"

The answer is: "Depends on the shape." So start discarding the smallest magnitude coefficients and see how much the accuracy suffers.

When your are constructing the Fourier series, you need to divide the unnormalized DFT coefficients by $N$. You also want to count the upper half of the DFT as negative frequency, so $N-1$ corresponds to $-1$, etc.

So basically you are taking the DFT of a discrete sequence, and then reconstructing an interpolation using the coefficients.

Hope this helps.

Ced


I'm putting this in my answer as I don't want to trigger shunting this conversation off to a chat space (a policy I disagree with BTW).

The issue at hand isn't simply can a closed figure be parameterized, the question is firmly set as an application of the summation of epicycles (you know, how planetary motion used to be modelled before Copernicus's change of reference frame). Yes, there are other ways to parameterize circular motion rather than sine and cosine, but they are clumsy.

There are also other ways to step back and parameterize the figure as a whole and there is no requirement that it be periodic. A Legendre basis comes to mind. It just so happens with the DFT approach it is inherently periodic.

In my opinion the OP thought it was cool (as do I) that you can draw an arbitrary figure (within limits) and was trying to understand how the concept of epicycles relates to the DFT.

Let's do a little math to make it clearer. Using conventional normalization and notation, the DFT is:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i n k \frac{2\pi}{N} } $$

Since the $x[n]$ are known, the $X[k]$ are also now known. Now, let's look at the inverse:

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i n k \frac{2\pi}{N} } $$

If we simply allow $n$ to be real valued and treat inverse DFT definition as a continuous equation we run into trouble past the Nyquist frequency. In the discrete case there is no difference because they will match at the sample points. In between, it does. So the summation has to be shifted to be centered around the DC bin. (Assuming $N$ is even, otherwise it can be similarly worked out)

$$ x[n] = \frac{1}{N} \sum_{k=-N/2+1}^{N/2} X[k] e^{i n k \frac{2\pi}{N} } $$

The equation can also be split into its real and imaginary parts:

$$ \Re(x(n)) = \sum_{k=-N/2+1}^{N/2} X[k] \frac{1}{N} \cos( n k \frac{2\pi}{N} ) $$ $$ \Im(x(n)) = \sum_{k=-N/2+1}^{N/2} X[k] \frac{1}{N} \sin( n k \frac{2\pi}{N} ) $$

I would have used $x$ and $y$, but $x$ is already taken. These equations are clearly in the form of Fourier series with coefficients of $ X[k] / N $.

I'm not trying to educate r b-j here, I know he knows this stuff thoroughly. I'm simply saying bringing in alternative parameterizations, or alternative coordinate systems, is a distraction from the core issues at hand.


Yep, just educated by r b-j. Thanks for the edits too.

Indeed the Nyquist term should be split in half and the result is those two epicycles will cancel each other's imaginary parts and double up the real part. Since there aren't a lot of zig zags in the figure, I would expect the magnitude of this coefficient to be low.

Just for kicks, I wrote a little Gambas program to demonstrate the math. You can find it here:

https://forum.gambas.one/viewtopic.php?f=4&t=725

I also did a little bit of a freehand Pi symbol. Sure, it looks a little drunk, but it still demonstrates the point.

enter image description here


Per request, here is a little bit of corner treatment. The corners worked better than I expected. I think this example truly exemplifies what I said earlier about the really interesting problem being finding the point placements along the figure that yield the closest fit.

enter image description here


What I haven't said explicitly in this discussion is that the complex value of $ e^{i\theta} $ moves along the complex unit circle, and is thus a model of an epicycle, so each of the products inside the loops represents the radius location of its respective epicycle (i.e. a line segment) at that point in time if you want to do the epicycle animation. The length of the radius is the magnitude of the coefficient as the magnitude of $ e^{i\theta} $ is always 1.

Complex.Polar( r, theta ) = r * e ^ { i theta }

You might find this article of mine helpful in understanding this material:

I'm not a fan of MATLAB (mostly because of the extremely nearsighted use of one based arrays), so I will refrain commenting on your pseudo-code. Instead, here is my code that actually calculates the value of the interpolation at a given "n".

You can follow the link and download it yourself (I just put the new version up that allows multiple figures in the same drawing). If you have Linux, you can install Gambas (PPA:gambas-team/gambas3) to run it.

[Note: the n in the code is just an iterator, the t is the actual n, I'm not bothering to edit the code.]

.
.
.

        For n = 0 To myPoints.Count * 100 - 1
          t = n / 100  
          p = Calculate(t, w)
          Paint.Arc(p.Real, p.Imag, 1) 
          Paint.Fill() 
        Next
.
.
.

'=======================================================================
Public Sub Calculate(ArgT As Float, ArgDFT As Vector) As Complex

        Dim k, N As Integer
        Dim p As Complex
        Dim a, b As Float

        N = ArgDFT.Count

        b = ArgT * Pi(2) / N

        If Even(N) Then
           GoSub EvenCase 
        Else    
           GoSub OddCase 
        Endif

        Return p

'-----------------------------------------------------------------------
EvenCase:

        p = ArgDFT[0] + ArgDFT[N / 2] * Cos(ArgT * Pi)

        For k = 1 To N / 2 - 1
          a = b * k
          p += ArgDFT[k] * Complex.Polar(1.0, a)
          p += ArgDFT[N - k] * Complex.Polar(1.0, -a)
        Next

        Return

'-----------------------------------------------------------------------
OddCase:

        p = ArgDFT[0]

        For k = 1 To (N - 1) / 2
          a = b * k
          p += ArgDFT[k] * Complex.Polar(1.0, a)
          p += ArgDFT[N - k] * Complex.Polar(1.0, -a)
        Next

        Return

End
'=======================================================================

A rebuttal to r b-j:

Robert, I strongly disagree with several of your assertions.

1) The interpolated points (and the path they form by LineTo calls) will follow whatever order you feed the points

2) Using a 0 to $2\pi$ range for "t" confuses the issue when compared to my answer in which "t" ranges from 0 to N, i.e. the same scale as the discrete scale, only including the real values in between the integers. This is a reflection of your frame of reference being the continuous case. [No longer relevant, I've ditched the t]

3) Treating (x,y) as a vector, rather than a single complex value x + iy, separates the parameterization into two independent problems which need not necessarily be parameterized by the same methodology. It is only in the complex value interpretation that the concept of Epicycles, which is the core of this problem, is meaningful.

4) Bunching the points in the corners, without a sufficient number of points in between, will cause overruns on the corners. The demonstration of this is why I added the fourth figure in my last graph.

5) Your definition of $a_k$ and $b_k$ is meaningless as there is no continuous function given, only a set a sample points. Therefore the Fourier coefficients should be calculated using the discrete definition, i.e. a summation not an integration. You have put the cart before the horse. With a different parameterization, like Legendre, you won't have a repeat pattern outside the range, won't necessarily match in between the points, but you will match at all the sample points.

The title question is: "How to get Fourier coefficients to draw any shape using DFT?"

The answer is: "The normalized DFT bin values are the Fourier coefficients."

In other words, simply replacing an integer "n" in the centered inverse DFT with a continuous real valued variable will produce the interpolated results. You can't get any more elegant than that. My code is an expression of this. I am assuming that the OP will implement it in MATLAB (with the necessary index adjustment).

You are making this way more complicated than it needs to be.


Here is an equivalent coding of the even case loop to clarify the meaning of "k" and it's range.

       For k = -N / 2 + 1 To -1
          a = b * k
          p += ArgDFT[k+N] * Complex.Polar(1.0, a)
       Next

       For k = 1 To N / 2 - 1
          a = b * k
          p += ArgDFT[k] * Complex.Polar(1.0, a)
       Next


This one is for Olli, using N = 9. If the figure is indeed a triangle, you can see with the proper point placement, a better fit can be found that does hit all the points as well. Of course, more points (more epicycles) could be added to get an even closer fit.

enter image description here

The auxiliary problem here (mentioned before), and I thought it would be the one you would tackle, Olli, is how to place the sample points on the underlying continuous figure to minimize either the "ripples" or the "overruns".


This is what happens when you treate the upper half of the DFT as positive frequencies rather than negative ones. You can clearly see that all the points get hit, but in between the results are not what is desired. Maybe there are some novelty applications where this would be beneficial.

enter image description here

I did this in response to Olli's challenge of if it could be done with positive frequencies only. Perhaps, if the real and imaginary are separated and cosine series are used for the two parameterizations, but I think that thwarts the intent of the question, and it wouldn't be an epicycle implementation anymore.

My initial inclination is to say no. I think the question is equivalent to "Can you construct a counter-clockwise corkscrew out of a summation of clockwise corkscrews?" Maybe with an infinite number, I've seen too many weird things in math concerning approaching infinity to rule it out, but here I can't even imagine a sequence that is an approximation.


Chris,

I have nothing against capital letters for variables. Indeed, I like to use $S_n$ for my signal values. What I don't like is using a lower case "x" for the signal and an upper case "X" for the DFT. To me, that isn't a sufficient distinction as they are describing two totally different domains. In addition "X" is one of letters that the lower case version and upper case version is most similar, making them even harder to distinguish in hand written math.

We basically have three scales (or function domains) in this situation:

1) n goes from 0 to N-1 on the integers for the sample points (for the input points and the output of the inverse DFT)

2) k goes from 0 to N-1 on the integers in the inverse DFT definition, then shifted half a frame to de-alias the upper half

3) t goes from 0 to $2\pi$ is the domain for the series solution (you and Robert) and (0 to N - 1/100 in my code)

So yeah, you are being misleading by using K in the T domain.

In my code, ArgDFT is the 1/N normalized DFT, and ArgT is my original "t" parameter, which is on the same scale as "n", but continuous. My "b" in the code is the same as your "t".

In summary of the process:

When you take the 1/N normalized DFT of a sample sequence you are simultaneously finding the coefficients for the continuous Fourier series which will pass through all the points. (A strong argument for why 1/N normalization should be the convention to use).

The domain of the series solution can be rescaled by a variable substitution:

$$ n = t \cdot \frac{N}{2\pi} $$

into the inverse DFT interpreted as a continuous function.

$$ x(n) = \frac{1}{N} \sum_{k} X[k] e^{i n k \frac{2\pi}{N} } $$

$$ z(t) = x(t \cdot \frac{N}{2\pi}) = \frac{1}{N} \sum_{k} X[k] e^{i t \cdot \frac{N}{2\pi} k \frac{2\pi}{N} } $$

$$ z(t) = \sum_{k} \frac{X[k]}{N} e^{i k t } $$

That is the series solution of the continuous interpolating path. It is just a function of t. You can differentiate it in respect to t to figure out your "pen velocity" if you want.

It is clear that you are now understanding what I meant by "the point placement problem", and it looks like Olli's interest has been piqued in it as well.

If you haven't already, I would suggest that you reread everything in this thread. With a better basis of understanding, the things that have been said should be more meaningful.


Epilogue: A different perspective, familiar to many here, on the situation. However, it doesn't yield the Fourier coefficients.

Meant to be a slog pile.

$$ z(t) = \sum_{k} \frac{X[k]}{N} e^{i k t } $$

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i n k \frac{2\pi}{N} } $$

$$ z(t) = \frac{1}{N} \sum_{k} \sum_{n=0}^{N-1} x[n] e^{-i n k \frac{2\pi}{N} } e^{i k t } $$

$$ z(t) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] \sum_{k} e^{i k ( t - \frac{n}{N}2\pi )} $$

$$ t_n = t - \frac{n}{N}2\pi $$

$$ D(t_n) = \sum_{k} e^{i k t_n } $$

$$ z(t) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] D(t_n) $$

Odd case: $k = -(N-1)/2 \to (N-1)/2$

Let $l = k + (N-1)/2$ goes $0 \to N-1$

$$ k = l - (N-1)/2 $$

$$ \begin{aligned} D(t_n) &= \sum_{l=0}^{N-1} e^{i ( l - (N-1)/2 ) t_n } \\ &= \sum_{l=0}^{N-1} e^{i l t_n } e^{-i \frac{N-1}{2} t_n } \\ &= e^{-i \frac{N-1}{2} t_n} \sum_{l=0}^{N-1} (e^{i t_n })^l \\ &= e^{-i \frac{N-1}{2} t_n} \frac{1 - (e^{i t_n })^N }{ 1 - e^{i t_n } } \\ &= e^{-i \frac{N-1}{2} t_n} \left[ \frac{e^{i t_n N / 2 } } { e^{i t_n / 2 } } \cdot \frac{ e^{-i t_n N / 2 } - e^{i t_n N/2 } }{ e^{-i t_n / 2 } - e^{i t_n / 2 } } \right] \\ &= \frac{e^{i t_n N / 2 } - e^{-i t_n N / 2 }} { e^{i t_n / 2 } - e^{-i t_n / 2 } } \\ &= \frac{ 2i \cdot \sin( N t_n / 2 ) } { 2i \cdot \sin( t_n / 2 ) } \\ &= \frac{ \sin( N t_n / 2 ) } { \sin( t_n / 2 ) } \end{aligned} $$

$$ z(t) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] \frac{ \sin( N t_n / 2 ) } { \sin( t_n / 2 ) } $$

$$ z(t) = \sum_{n=0}^{N-1} x[n] \frac{ \sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \sin( (t - \frac{n}{N}2\pi) / 2 ) } $$

Notice that the quotient is real valued so it can be thought of as a weight and the summation is then the time variant weighted average of the set of sample points.


Epilogue II

After a whole lot of discussion on other questions, it is a apparent that the Nyquist bin should be split evenly between the negative and positive frequencies.

Even case: $k = 1/2 ( N/2 \text{ and } -N/2 ), -N/2 + 1 \to N/2 - 1 $

Let $l = k + N/2 - 1 $ goes $0 \to N-2$

$$ k = l - N/2 + 1 $$

$$ \begin{aligned} D(t_n) &= \frac{1}{2} \left[ e^{i ( N/2 ) t_n } + e^{i (-N/2 ) t_n } \right] + \sum_{l=0}^{N-2} e^{i ( l - N/2 + 1 ) t_n } \\ &= \cos \left( \frac{N}{2} t_n \right) + \sum_{l=0}^{N-2} e^{i l t_n } e^{i (- N/2 + 1 ) t_n } \\ &= \cos \left( \frac{N}{2} t_n \right) + e^{i (- N/2 + 1 ) t_n } \sum_{l=0}^{N-2} (e^{i t_n })^l \\ &= \cos \left( \frac{N}{2} t_n \right) + e^{i (- N/2 + 1 ) t_n } \frac{1 - (e^{i t_n })^{N-1} }{ 1 - e^{i t_n } } \\ &= \cos \left( \frac{N}{2} t_n \right) + e^{i (- N/2 + 1 ) t_n } \left[ \frac{e^{i t_n ( N - 1 ) / 2 } } { e^{i t_n / 2 } } \cdot \frac{ e^{-i t_n ( N - 1 ) / 2 } - e^{i t_n ( N - 1 ) / 2 } }{ e^{-i t_n / 2 } - e^{i t_n / 2 } } \right] \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{e^{i t_n ( N - 1 ) / 2 } - e^{-i t_n ( N - 1 ) / 2 }} { e^{i t_n / 2 } - e^{-i t_n / 2 } } \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{ 2i \cdot \sin( t_n ( N - 1 ) / 2 ) } { 2i \cdot \sin( t_n / 2 ) } \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{ \sin( t_n N /2 ) \cos( t_n / 2 ) - \cos( t_n N /2 ) \sin( t_n / 2 ) } { \sin( t_n / 2 ) } \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{ \sin( t_n N /2 ) } { \sin( t_n / 2 ) } \cos( t_n / 2 ) - \cos( t_n N /2 ) \\ &= \frac{ \sin( N t_n/2 ) }{ \sin( t_n / 2 ) } \cos( t_n / 2 ) \end{aligned} $$

The above derivation can be done using coefficients other than 1/2 and 1/2 for the positive and negative Nyquist terms, but then the simplification that occurs towards the end wouldn't happen and the expression would be more complicated. It would also be abundantly clear that if the set of $x[n]$ were real the interpolation wouldn't necessarily be real. For 1/2 and 1/2, the interpolation values will all be real.

The continuous interpolation function is then:

$$ \begin{aligned} z(t) &= \frac{1}{N} \sum_{n=0}^{N-1} x[n] \left[ \frac{ \sin( N t_n / 2 ) }{ \sin( t_n / 2 ) } \right] \cos( t_n / 2 ) \\ &= \sum_{n=0}^{N-1} x[n] \left[ \frac{ \sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \sin( (t - \frac{n}{N}2\pi) / 2 ) } \right] \cos( (t - \frac{n}{N}2\pi) / 2 ) \\ &= \sum_{n=0}^{N-1} x[n] \frac{ \sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \tan( (t - \frac{n}{N}2\pi) / 2 ) } \end{aligned} $$

It is quite remarkable that this formula matches the odd case version with a simple "window function" applied as seen in the first two lines. The last matches R B-J's given formula which is in more of a concise format.

Looking at the case where N = 2

$$ \begin{aligned} z(t) &= x[0] \left[ \cos^2( t / 2 ) \right] + x[1] \left[ \cos^2( (t - \pi) / 2 ) \right] \\ &= x[0] \left[ \frac{ \cos( t ) + 1 }{2} \right] + x[1] \left[ \frac{ \cos( t - \pi ) + 1 }{2} \right] \\ &= \frac{1}{2} ( x[0] + x[1] ) + \frac{1}{2} ( x[0] - x[1] ) \cos( t ) \end{aligned} $$

Which means the alternating sequence 1, -1, 1, -1, gets interpolated as:

$$ z(t) = \cos( t ) $$

Which is a direct result of the Nyquist split assumption.

$\endgroup$
  • $\begingroup$ Ced, are they my comments you don't understand? $\endgroup$ – robert bristow-johnson Jun 24 at 9:37
  • $\begingroup$ OK I could try to DFT with my existing unequally spaced samples but I've got an electrical engineering background and my (limited) understanding is tightly tied to frequency bins equally spaced by $F_s/N$ where $F_s$ is sampling frequency. So I have to ask - is there a way to convert my unequally samples to equally spaced? For example, interpolating between existing samples and sampling that new function at equal spacing. $\endgroup$ – Chris Jun 24 at 9:44
  • $\begingroup$ Actually tonight I made some promising progress along the suggested lines using my unequally spaced samples, I just need some time to work it through. Meanwhile thanks ced and robert for clues so far, I'll be back soon ... $\endgroup$ – Chris Jun 24 at 11:54
  • $\begingroup$ @robertbristow-johnson Yes, I don't understand why you brought polar representation or parametric equations into it. This is merely an exercise of using a DFT for interpolation (demonstrating that an inverse DFT can be interpreted as a sum of epicycles, e.g. Euler's equation). What the OP may not realize yet is that it will not match his figure between the points exactly, but will pass through every one of the sample points. $\endgroup$ – Cedron Dawg Jun 24 at 12:32
  • 3
    $\begingroup$ there is such a facility. i really dislike it, but sometimes the mods come around and "encourage" us to take it to chat. @PeterK, you ain't gonna do that, are you? $\endgroup$ – robert bristow-johnson Jun 27 at 0:00
2
$\begingroup$

Complex Fourier series of a piece-wise linear waveform tracing the desired shape

Instead of using discrete Fourier transform (DFT) / fast Fourier transform (FFT), a more direct approach is to define a piece-wise linear continuous-time waveform that traces the desired shape on the complex plane, and to directly calculate its Fourier series. Bezier curves or such could be used for shape definition and approximated using line segments to arbitrary accuracy. Your third figure is already sketched using line segments. We can use its node (corner) coordinates, but the times of the nodes of the waveform will need to be pulled out of a hat. We shall go with uniform sampling in time, without hard-coding that in the math or the scripts. The waveform can be plotted in Octave, taking as input the node coordinates from the file pi.csv from this answer:

graphics_toolkit("gnuplot")  # Octave specific to get prettier plots

xy = csvread("pi.csv");
z = xy(:,1) + i*xy(:,2);
M = length(xy);
t = (0:M-1)'*2*pi/M;

plot([t; 2*pi], [real(z); real(z(1))], "b");
hold on
plot([t; 2*pi], [imag(z); imag(z(1))], "r");
plot(t, real(z), "k.");
plot(t, imag(z), "k.");
xlim([0,2*pi])
ylim([-250,250])
xlabel("t")
hold off

![enter image description here
Figure 1. The real (blue) and imaginary (red) parts of a piece-wise linear waveform tracing the desired shape.

Complex Fourier series of the waveform

Let's have a look at a single linear segment of the waveform. A $2\pi$-periodic continuous-time waveform that is otherwise zero but has a line segment starting at complex value $z_0 = x_0 + y_0\,i$ at time $t_0$ and ending with value $z_1 = x_1 + y_1\,i$ at time $t_1 > t_0$ has coefficients of its complex Fourier series (using your second equation):

$$c_k=\frac{1}{2\pi}\int_\limits{t_0}^{t_1}\left(z_0 + \frac{t-t_0}{t_1-t_0}(z_1 - z_0)\right)\, e^{-ikt} \mathop{dt},\tag{1}$$

where $\frac{t-t_0}{t_1-t_0}$ goes from $0$ to $1$ as $t$ goes from $t_0$ to $t_1$. For $k = 0$ we have:

$$c_0 = \frac{(t_1 - t_0)(z_0 + z_1)}{4\pi},\tag{2}$$

and for negative and positive $k \ne 0$:

$$\begin{gather}c_k = \color{blue}{\frac{z_1\,\sin(k\,t_1) - z_0\,\sin(k\,t_0)}{2\pi k}} + \frac{(z_1 - z_0) \cos(k\,t_1) - (z_1 - z_0) \cos(k \,t_0)}{2\pi k^2\,(t_1 - t_0)}\,+\\ i\,\left(\color{blue}{\frac{z_1\,\cos(k\,t_1) - z_0\,\cos(k\,t_0)}{2\pi k}} - \frac{(z_1 - z_0)\,\sin(k\,t_1) - (z_1 - z_0) \sin(k\,t_0)}{2\pi k^2\,(t_1 - t_0)}\right).\end{gather}\tag{3}$$

You'd fully populate the range $0 \le t \le 2\pi$ (or any range of length $2\pi$, for example $-\pi \le t \le \pi$ compatible with your second equation) with non-overlapping linear segments, and for each integer $k$ separately, use the sum of the coefficients over the segments to obtain the corresponding coefficient for the complex Fourier series of the full piece-wise linear waveform. This works, because addition in frequency domain is equivalent to addition in time domain which splices the linear segments together. The terms colored in blue in Eq. 3 will cancel in the sum over waveform segments and need not be included in it. Using your first equation, the full waveform equals its complex Fourier series:

$$z(t) = \sum_{k=-N}^N c_k\,e^{ikt},\tag{4}$$

with $N=\infty$. Except for the constant term with coefficient $c_0$, the individual harmonic terms of the sum go in circles on the complex plane, and can be considered epicycles.

You can truncate the series at some finite $N$. Here is an Octave script implementing this approach. It calculates a truncated complex Fourier series of the piece-wise linear waveform defined by the given nodes (corners), assuming uniform time distribution of nodes:

graphics_toolkit("gnuplot")  # Octave specific to get prettier plots

xy = csvread("pi.csv");
z = xy(:,1) + i*xy(:,2);

M = length(xy);
N = floor(M/2) - 1;  # Truncation length, this can be any positive integer
k = -N:N;
t = (0:M-1)'*2*pi/M;  # This can be any ascending sequence of times of the nodes obeying 0 <= t < 2 pi
t1 = circshift(t,-1);
t1(end) = 2*pi;
z1 = circshift(z,-1);

c = sum(((z1 - z).*cos(k.*t1) - (z1 - z).*cos(k .*t))./(2*pi*k.^2.*(t1 - t))+ i*(-((z1 - z).*sin(k.*t1) - (z1 - z).*sin(k.*t))./(2*pi*k.^2.*(t1 - t))), 1);
c(N + 1) = sum(((t1 - t).*(z + z1))/(4*pi), 1);

# c now contains complex Fourier series coefficients in order k

z_new = (2*N + 1)*ifft(ifftshift(c));  # Uniformly sample the reconstruction in time
xy_new = [real(z_new)', imag(z_new)'];
csvwrite("pi_bandlimited.csv", xy_new);  # Save samples. This should work with rb-j's script for odd length

os = 8;  # Oversampling factor, integer
z_os = os*(2*N + 1)*ifft([c(N+1:end) zeros(1, (N*2+1)*(os - 1)) c(1:N)]);  # Band-limited approximation
plot([real(z_os) real(z_os(1))], [imag(z_os) imag(z_os(1))], "-")
xlim([-250,250]);
ylim([-250,250]);
hold on
#plot(real(z_new), imag(z_new), "+")  # New samples
#plot([real(z);real(z(1))], [imag(z);imag(z(1))], "-")  # Desired shape
hold off

The result (Fig. 2) can be visually compared against that from a Fourier interpolation (given by a DFT-based approach presented in @robertbristow-johnson's answer) (Fig. 3), continuing the above Octave script:

z_ftos = interpft(z, length(z)*os);  # Fourier interpolate
plot([real(z_ftos); real(z_ftos(1))], [imag(z_ftos); imag(z_ftos(1))], "-")
xlim([-250,250]);
ylim([-250,250]);

enter image description here
Figure 2. The result of the approach suggested in this answer is a least-squares band-limited approximation of the piece-wise linear waveform tracing desired shape, here using 119 harmonic terms.

enter image description here
Figure 3. The result of the Fourier interpolation approach (not presented in this answer), using 120 harmonic terms.

As can be seen in the above, the suggested approach results in cleaner tracing of the desired shape, and can be easily made more accurate by increasing $N$ (Fig. 4).

enter image description here
Figure 4. The result of the suggested approach setting N = 2000 truncates the complex Fourier series to 4001 harmonic terms and gives this visually indistinguishable approximation of the desired shape.

If you wish to, you can rearrange the sum of Eq. 4 to interleaved positive and negative coefficients, or order it by decreasing $|c_k|$. You can also construct a sparse approximation by picking only the largest coefficients. We can have a look at the magnitudes of the harmonics by running the above Octave script with N = 20000 and plotting:

loglog(abs(k), abs(c), '.');
xlim([1, 20000]);
ylim([0.000001, max(abs(c))]);
xlabel("|k|");
ylabel("|c_k|")

enter image description here
Figure 5. Magnitudes of the harmonics of the piece-wise linear waveform tracing the $\pi$ shape. Due to the continuity of the waveform, the envelope decays asymptotically with a -40 dB / decade slope.

Waveform approximation error

The mean square error in the approximation of the waveform will go down with every included term of the complex Fourier series, no matter which order. This is because the harmonic terms are orthogonal so the mean square of any partial sum of the terms is a sum of the mean squares of the terms, which are:

$$\frac{1}{2\pi}\int_\limits{0}^{2\pi}\left|c_k\,e^{ikt}\right|^2\mathop{dt} = |c_k|^2.\tag{5}$$

The mean square of the waveform equals the mean square of the complex Fourier series which is the limit of the partial sum of the harmonic terms as $N\to\infty$ and can be equivalently calculated in time domain as the sum of the mean squares of the linear segments, which are:

$$\frac{1}{2\pi}\int_\limits{t_0}^{t_1}\left|z_0 + \frac{t-t_0}{t_1-t_0}(z_1 - z_0)\right|^2 \mathop{dt} = \frac{(t_1 - t_0)(z_0^2 + z_0 z_1 + z_1^2)}{6\pi}.\tag{6}$$

The difference between the two sums is the mean square error of the waveform approximation, which however is not a uniquely defined or even always a reasonable metric of error in tracing the desired shape.

Optimal tracing of the shape

A truncated complex Fourier series of a piece-wise linear complex waveform is usually not the optimal $2N+1$-term approximation that most closely traces the desired shape. It is possible to construct a new set of line segments which is a time-stretched version of the complex waveform that will trace the same shape but can result in a better truncated series. I think it can be made arbitrarily close to in some sense the optimal waveform for that shape if the number of line segments is increased sufficiently. Finding the optimal waveform seems like a difficult problem to tackle.

$\endgroup$
  • $\begingroup$ No, Olli, your are missing the point(s), literally. All you have done is produce an approximation of linear interpolation that is epicycle implementable, taking roughly twice as many epicycles to do so. I will give a fuller explanation later about the fundamental flaw in both yours and Robert's approaches. $\endgroup$ – Cedron Dawg Jun 29 at 10:34
  • $\begingroup$ Strike the twice as many portion of that last comment, you are taking even more than that. The limit as N goes to infinity is the linear interpolation which I am pretty sure is not anybody's idea of the best solution. $\endgroup$ – Cedron Dawg Jun 29 at 10:42
  • $\begingroup$ Yes, but your method will not include the sample points except for special cases until an infinite number of terms are taken. I think your characterization of it being "a more direct approach" is completely off base. From a pragmatic programming aspect, your method is a nightmare in comparison. Although not explicitly specified, I do think the OP was looking for a solution that would give a nice curve through the points rather than line segments. However, in the title question itself, it was explicitly indicated that the OP wanted to know how to get the coefficients from the DFT. $\endgroup$ – Cedron Dawg Jun 29 at 11:47
  • $\begingroup$ I believe the OP is using linear interpolation because it is the simplest method to "connect the dots" to draw a figure. Did you see the triangle samples I added at the bottom of my post for you? $\endgroup$ – Cedron Dawg Jun 29 at 12:13
  • $\begingroup$ And the challenge? $\endgroup$ – Cedron Dawg Jun 29 at 12:14
2
$\begingroup$

Input data files

This answer is for storing input data files for testing solutions to the problem.

Here are the 120 $x,y$ coordinates for the $\pi$ shape, by @Chris. Save as pi.csv:

108,0
110,25
112,50
113.5,75
115,100
116,125
117.5,150
125,150
150,150
175,150
200,150
225,150
225,175
225,200
225,220
200,220
175,220
150,220
125,220
100,220
75,220
50,220
25,220
0,219.5
-25,219
-50,217
-75,215
-100,212
-125,209
-150,203
-158,200
-175,190
-190,175
-203,150
-211,125
-220,100
-225,85
-209,85
-200,100
-182,125
-175,132
-150,145
-125,150
-100,150
-87,150
-87.5,125
-89,100
-92,75
-95,50
-100,25
-105,0
-113,-25
-122,-50
-136,-75
-152,-100
-170,-125
-186,-150
-189,-175
-178,-200
-175,-205
-150,-220
-125,-220
-100,-202
-85,-175
-77,-150
-73,-125
-70,-100
-67.5,-75
-65,-50
-62,-25
-60,0
-57,25
-54.5,50
-51.5,75
-49,100
-47,125
-45,150
-25,150
0,150
25,150
50,150
58,150
55,125
53,100
51,75
49,50
47,25
44.5,0
42,-25
40,-50
38.5,-75
37.5,-100
37,-125
37.5,-150
43,-175
49,-185
66,-200
75,-205
100,-215
125,-218
150,-214
175,-203
179,-200
201.5,-175
213,-150
221,-125
226.5,-100
227.5,-88
210,-88
209,-100
200,-123
197,-125
175,-141
150,-144
125,-134
117,-125
109,-100
106,-75
106,-50
106.5,-25
$\endgroup$
  • $\begingroup$ By Jove, I think he's got it. Nicely done, Chris. The next challenge is to achieve your (or any) figure with the fewest possible epicycles. $\endgroup$ – Cedron Dawg Jul 1 at 19:00
  • $\begingroup$ @Cedron ... um .... this was Olli's work, not mine ... yet. But I'll take the kudos for the hack work behind providing the sample points :) I'm still slogging my way through all the comments and formulas from days ago. $\endgroup$ – Chris Jul 1 at 23:23
  • $\begingroup$ @Chris, Oops, missed that. Olli always does spectacular work which is generally quite thorough. Hopefully your slogging isn't too tough. $\endgroup$ – Cedron Dawg Jul 2 at 14:05
  • $\begingroup$ @Chris I moved your modified r b-j's scripts and the result plots to the question, because the modifications broke the original functionality. The plots seem to only show data loaded from files. pi.csv as defined here has 120 rows as it defines a closed curve without explicitly duplicating the first point at the end as a 121th row, which is something that can be done in the script that loads it, if needed. I restored r b-j's original scripts into his answer. $\endgroup$ – Olli Niemitalo Jul 2 at 14:44
  • $\begingroup$ @Olli I have responded in the chat session for r b-j's answer. I deliberately removed the 121st point as it is simply a repeat of the first. point. If that's wrong let's talk about it in the chat session. $\endgroup$ – Chris Jul 2 at 22:53
0
$\begingroup$

I know that I am not really being listened to, Chris, but I know exactly what you're trying to do. I know exactly what the problem is. I know exactly what the mathematics are. And I know exactly what you should do and exactly how you should think about it. You're starting to move in the right direction with an ordered set of $N$ points with a horizontal component $x_n$ and a vertical component, $y_n$. And you are expressing these pairs together as a set of complex numbers:

$$z_n \triangleq x_n + i y_n \qquad \qquad \text{for } 0 \le n < N $$

For the moment, let's toss MATLAB (and return to it later), because of their horrible off-by-one problem of indexing. That can be fixed later simply by recognizing the indexes are off by one.

Now, because this is a curve in only two-dimensions, you can represent it as a curve on the complex plane $\mathbb{C}$, but you really should be thinking of the points as simple coordinate pairs of real numbers $(x_n, y_n)$ in the two-dimensional real space $\mathbb{R}^2$ and get two Fourier series for the two real functions $x(t)$ and $y(t)$, since these two real functions are periodic having the same period and they are synchronized together by definition. When $x(t) = x_n$ that will happen at the very same "time" $t$ that $y(t)=y_n$. This is why the points need not be equally spaced. But they do need to be in order, with no skipping over any points.

You can do this with a single set of complex values $z_n$ in $\mathbb{C}$ and it will lose you nothing nor gain you nothing except maybe you can perform just one FFT instead of two.

Start anywhere on the curve and call that point "$(x_0, y_0)$" or the complex "$z_0 = x_0 + i y_0$" if you must. Then move (I would suggest) in the counter-clockwise direction and the very next point you hit should be $(x_1, y_1)$ or $z_1$ and continuing counter-clockwise the next point is $(x_2, y_2)$ or $z_2$. Do not skip over points. Eventually you will go around the entire simple closed curve and get to the last point $(x_{(N-1)}, y_{(N-1)})$ or $z_{(N-1)}$ and then the "pen" moves further and you get right back to your starting location of $(x_0, y_0)$ or $z_0$. And if you want this to look good, you will put a higher density of points around the sharp corners than elsewhere on the curve with lower curvature.

What you are doing is bandlimited reconstruction of two periodic waveforms, $x(t)$ and $y(t)$, both with periods of $2\pi$ from the two periodic sequences $x_n$ and $y_n$, both having a common period of $N$, that is you can think of your finite length sequences of $x_n$ and $y_n$ as one period of an infinitely long periodic sequence:

$$ x_{(n+N)} = x_n \qquad \qquad \forall n \in \mathbb{Z} $$

$$ y_{(n+N)} = y_n \qquad \qquad \forall n \in \mathbb{Z} $$

Likewise the two periodic functions $x(t)$ and $y(t)$ satisfy:

$$ x(t+2\pi) = x(t) \qquad \qquad \forall t \in \mathbb{R} $$

$$ y(t+2\pi) = y(t) \qquad \qquad \forall t \in \mathbb{R} $$

and from Fourier, we know that,

$$ x(t) = \sum\limits_{k=-\infty}^{\infty} a_k \, e^{i k t} $$

$$ y(t) = \sum\limits_{k=-\infty}^{\infty} b_k \, e^{i k t} $$

where

$$ a_k = \frac{1}{2\pi} \int\limits_{t_0}^{t_0+2\pi} x(t) \, e^{-ikt} \, \mathrm{d}t \qquad \forall t_0 \in \mathbb{R} $$

and

$$ b_k = \frac{1}{2\pi} \int\limits_{t_0}^{t_0+2\pi} y(t) \, e^{-ikt} \, \mathrm{d}t \qquad \forall t_0 \in \mathbb{R} $$

If $x(t)$ and $y(t)$ are purely real, then the Fourier coefficients with negative indices are complex conjugate of their counterparts having positive indices:

$$ a_{-k} = a_k^* $$ $$ b_{-k} = b_k^* $$ This is called "Hermitian symmetry". _

Of course you can define:

$$ z(t) \, \triangleq \, x(t) \, + \, i \, y(t) $$

and say that:

$$ z(t+2\pi) = z(t) \qquad \qquad \forall t \in \mathbb{R} $$

$$ z(t) = \sum\limits_{k=-\infty}^{\infty} c_k \, e^{i k t} $$

$$ c_k = \frac{1}{2\pi} \int\limits_{t_0}^{t_0+2\pi} z(t) \, e^{-ikt} \, \mathrm{d}t \qquad \forall t_0 \in \mathbb{R} $$

where the Fourier coefficients $c_k = a_k \, + \, i \, b_k$ and it may seem elegant. Because $z(t)$ is not purely real, there do not have Hermitiam symmetry and we cannot say that $c_{-k}$ is the complex conjugate of $c_k$.


Now the Discrete Fourier Transform takes one period of length $N$ of a periodic sequence, $x_n$ (DSPers like to say "$x[n]$" and we like "$j$" instead of "$i$") and transforms it into one period, also of length $N$, of another period sequence, we'll call $X_k$ (or "$X[k]$").

$$ X_k \triangleq \sum\limits_{n=0}^{N-1} x_n e^{-i 2 \pi \frac{nk}{N}} $$

the inverse DFT transforms it back:

$$ x_n = \frac{1}{N} \sum\limits_{k=0}^{N-1} X_k e^{+i 2 \pi \frac{nk}{N}} $$

Now because both $x_n$ and $X_k$ are periodic with period $N$;

$$ x_{(n+N)} = x_n \qquad \qquad \forall n \in \mathbb{Z} $$

$$ X_{(k+N)} = X_k \qquad \qquad \forall k \in \mathbb{Z} $$

then the DFT and iDFT can be expressed as:

$$ X_k \triangleq \sum\limits_{n=n_0}^{n_0+N-1} x_n e^{-i 2 \pi \frac{nk}{N}} \qquad \qquad \forall n_0 \in \mathbb{Z} $$

$$ x_n = \frac{1}{N} \sum\limits_{k=k_0}^{k_0+N-1} X_k e^{+i 2 \pi \frac{nk}{N}} \qquad \qquad \forall k_0 \in \mathbb{Z} $$

Textbooks don't usually say this, but it's true. But it's true only for integer $n$ or $k$.

So let's consider a real sequence $x_n$ (as it is above as the $x$-coordinate of each point). Then Hermitian symmetry applies and $a_{-k} = a_k^*$. So, first, let's consider $N$ to be odd. That means that $\frac{N}2$ does not exist as integer and that there is no "Nyquist component": $X_{N/2}$.

So let's look at the iDFT and set $k_0 = -\frac{N-1}{2}$:

$$ x_n = \frac{1}{N} \sum\limits_{k=-(N-1)/2}^{(N-1)/2} X_k e^{+i 2 \pi \frac{nk}{N}} $$

Now comparing to the continuous and periodic function $x(t)$ and scaling the period $2\pi$ to the discrete period $N$, we can set:

$$ x(t)\bigg|_{t=\frac{2\pi}{N}n} = x_n \qquad \qquad \forall n \in \mathbb{Z} $$

This could be thought of as essentially uniformly sampling $x(t)$ at equally spaced sammpling instances of $t=\frac{2\pi}{N}n$ and the sampling period (not the same as the "period of $x(t)$) is $\frac{2\pi}{N}$. $N$ times that sampling period is the period of $x(t)$. If we make that substution for $t$ in $x(t)$ we have:

$$\begin{align} x(t)\bigg|_{t=\frac{2\pi}{N}n} &= \sum\limits_{k=-\infty}^{\infty} a_k \, e^{i k t}\bigg|_{t=\frac{2\pi}{N}n} \\ &= \sum\limits_{k=-\infty}^{\infty} a_k \, e^{i 2\pi \frac{kn}{N}} \\ &= \sum\limits_{k=-(N-1)/2}^{(N-1)/2} \frac{1}{N} X_k \, e^{i 2\pi \frac{kn}{N}} \\ &= x_n \\ \end{align}$$

if we relate

$$ a_k = \begin{cases} \frac{1}{N} X_k \qquad & |k| \le \frac{N-1}{2} \\ 0 \qquad & |k| > \frac{N-1}{2} \\ \end{cases} $$

Now, since in the computer program, it's unlikely your DFT will have negative indices and, relying on periodicity, then this comes out as

$$ a_k = \begin{cases} \frac{1}{N} X_k \qquad & 0 \le k \le \frac{N-1}{2} \\ \frac{1}{N} X_{(k+N)} \qquad & -\frac{N-1}{2} \le k < 0 \\ 0 \qquad & |k| > \frac{N-1}{2} \\ \end{cases} $$

For even $N$, it's nearly the same, but the component that sits on the Nyquist frequency must be divided by 2 into a positive frequency component and a negative frequency component.

$$\begin{align} x(t)\bigg|_{t=\frac{2\pi}{N}n} &= \sum\limits_{k=-\infty}^{\infty} a_k \, e^{i k t}\bigg|_{t=\frac{2\pi}{N}n} \\ &= \sum\limits_{k=-\infty}^{\infty} a_k \, e^{i 2\pi \frac{kn}{N}} \\ &= \frac{1}{N} X_{N/2} e^{i \pi n} + \sum\limits_{k=-\frac{N}{2}+1}^{\frac{N}{2}-1} \frac{1}{N} X_k \, e^{i 2\pi \frac{kn}{N}} \\ &= \frac{1}{2N} X_{N/2} e^{-i \pi n} + \frac{1}{2N} X_{N/2} e^{i \pi n} + \sum\limits_{k=-\frac{N}{2}+1}^{\frac{N}{2}-1} \frac{1}{N} X_k \, e^{i 2\pi \frac{kn}{N}} \\ &= x_n \\ \end{align}$$

Note that the two terms with $e^{i \pi n}=(-1)^n$ whether it's "$i$" or "$-i$" in the exponent. This is the Nyquist component that is oscillating at exactly half of the sample rate. Now the Fourier series coefficients are:

$$ a_k = \begin{cases} \frac{1}{N} X_k \qquad & |k| \le \frac{N}{2}-1 \\ \frac{1}{2N} X_k \qquad & |k| = \frac{N}{2} \\ 0 \qquad & |k| > \frac{N}{2} \\ \end{cases} $$

or

$$ a_k = \begin{cases} \frac{1}{N} X_k \qquad & 0 \le k \le \frac{N}{2}-1 \\ \frac{1}{N} X_{(k+N)} \qquad & -\frac{N}{2}+1 \le k < 0 \\ \frac{1}{2N} X_{N/2} \qquad & |k| = \frac{N}{2} \\ 0 \qquad & |k| > \frac{N}{2} \\ \end{cases} $$

That is how you get the Fourier series coefficients for a bandlimited periodic function from the DFT. We can do exactly the same song-and-dance for the real $y(t)$ with samples $y_n$ Fourier series coefficients $b_k$ and DFT coefficients $Y_k$.

We can even do the same for $z(t) \triangleq x(t) + i y(t)$ with points

$$ z(t)\bigg|_{t=\frac{2\pi}{N}n} = z_n \qquad \qquad \forall n \in \mathbb{Z} $$

because the Fourier series (and that integral) are linear.

$$ z(t) = \sum\limits_{k=-\infty}^{\infty} c_k \, e^{i k t} $$

$$ c_k = \frac{1}{2\pi} \int\limits_{t_0}^{t_0+2\pi} z(t) \, e^{-ikt} \, \mathrm{d}t \qquad \forall t_0 \in \mathbb{R} $$

That means, simply, that the Fourier series coefficients are

$$ c_k \, = \, a_k \, + \, i \, b_k $$

It's nothing more sophisticated than that. Complex $z_n$ is okay, but you could just as well do it with points in 2-space $\big(x_n,y_n \big)$ and the mathematics will be no different except now there is no Hermitian symmetry. We cannot say that $c_k$ and $c_{-k}$ are complex conjugates of each other. And you only need to run a single DFT of $z_n$ instead of two separate DFTs on $x_n$ and $y_n$.

Other than that, the mathematics are exactly the same. Pick $\big(x_n,y_n \big)$ or pick $z_n$, I don't care.


Here is the MATLAB program (Octave compatible) that r b-j kludged together to draw the $\pi$ shape by @Chris using coordinates from pi.csv from this answer:

clear;

xy = csvread("pi.csv");

NN = length(xy);
N = 65536;

xx = xy(:,1);
yy = xy(:,2);

figure(1);
plot(xx, yy, 'mx');

XX = fft(xx)/NN;
YY = fft(yy)/NN;

X = [XX(1+0:1+(NN/2-1)); 0.5*XX(1+NN/2); zeros(N-NN-1, 1); 0.5*XX(1+NN/2); XX(1+(NN/2+1):1+(NN-1))];
Y = [YY(1+0:1+(NN/2-1)); 0.5*YY(1+NN/2); zeros(N-NN-1, 1); 0.5*YY(1+NN/2); YY(1+(NN/2+1):1+(NN-1))];

x = N*ifft(X);
y = N*ifft(Y);

figure(1);
hold on;
plot(x, y, 'b');
hold off;

here is the result:

enter image description here

here is the same, but with one point removed so that NN is odd. note that there is no Nyquist value to split into two.

clear;

xy = csvread("pi.csv");
xy = vertcat(xy(1:111,:), xy(113:end,:));

NN = length(xy);
N = 65536;

xx = xy(:,1);
yy = xy(:,2);

figure(1);
plot(xx, yy, 'mx');

XX = fft(xx)/NN;
YY = fft(yy)/NN;

%   X = [XX(1+0:1+(NN/2-1)); 0.5*XX(1+NN/2); zeros(N-NN-1, 1); 0.5*XX(1+NN/2); XX(1+(N/2+1):1+(N-1))];
%   Y = [YY(1+0:1+(NN/2-1)); 0.5*YY(1+NN/2); zeros(N-NN-1, 1); 0.5*YY(1+NN/2); YY(1+(NN/2+1):1+(NN-1))];

X = [XX(1+0:1+((NN-1)/2)); zeros(N-NN, 1); XX(1+((NN+1)/2):1+(NN-1))];
Y = [YY(1+0:1+((NN-1)/2)); zeros(N-NN, 1); YY(1+((NN+1)/2):1+(NN-1))];

x = N*ifft(X);
y = N*ifft(Y);

figure(1);
hold on;
plot(x, y, 'b');
hold off;
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.