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.
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.
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.
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.
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.
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 :)
