FFT gives wrong answer for a generating function

I am trying to use the generating function below to make a symmetric Toeplitz matrix (a matrix where each diagonal has a constant value). Each is an element in a size n vector which forms the first row of my matrix A.

where: from

I believe that the terms are the Fourier Coefficients of and that I should be able to use a Fourier transform to calculate them. Specifically, I want to use a Fast Fourier Transform with fft().

The problem is that when I use fft(), I get coefficients that are both very large (~10^3) and complex. Since is real and even in the given interval, the Fourier coefficients should only be real.

Here is my code for how I am trying to use the FFT:

n = 50;
f = linspace(-pi, pi, n);   % Establish domain
f = f.^4 + 1;               % Calculate f(theta)
f = fft(f);                 % FFT

The output of f at the end is:

f =

   1.0e+03 *

   1.1047 + 0.0000i
   0.8360 + 0.0526i
   0.4487 + 0.0567i
   etc...

When I do this analytically, I get expected results:

f =

   20.4818
  -15.4784
    8.3696
    etc...

My question is: what am I doing wrong with the fft() function to give me such weird results? How can I fix it to give me the actual Fourier coefficients?

[Jan]

How did you solve this analytically? The problem is hidden there.

According to the documentation, FFT calculates:

        N
X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
       n=1

Especially X(1) equals sum(f), so 1.1047e3 is the expected result.

[Paul]

@Paul Anhalt

If you don't mind my asking, do you have to use the FFT for this problem? The FFT only yields an approximation to the alpha(k), but the exact form of alpha(k) is readily computed as shown below.

Answer 1

[Star Strider]

It is possible to get close to the desired result by ‘adjusting’ the fft arguments —

n = 50;
f = linspace(-pi, pi, n);   % Establish domain
f = f.^4 + 1;               % Calculate f(theta)
f = fft(f,(n-1))/(n-1)      % FFT              NOTE: Need To Divide By 'n'

f = 
  20.5359 + 0.0000i  15.5325 + 0.0000i   8.4239 + 0.0000i   4.1447 + 0.0000i   2.4285 - 0.0000i   1.5960 - 0.0000i   1.1338 + 0.0000i   0.8520 + 0.0000i   0.6680 - 0.0000i   0.5416 - 0.0000i   0.4512 - 0.0000i   0.3845 + 0.0000i   0.3341 - 0.0000i   0.2953 - 0.0000i   0.2648 - 0.0000i   0.2408 - 0.0000i   0.2215 + 0.0000i   0.2061 + 0.0000i   0.1938 + 0.0000i   0.1840 + 0.0000i   0.1764 - 0.0000i   0.1706 - 0.0000i   0.1664 - 0.0000i   0.1637 - 0.0000i   0.1623 - 0.0000i   0.1623 + 0.0000i   0.1637 + 0.0000i   0.1664 + 0.0000i   0.1706 + 0.0000i   0.1764 + 0.0000i   0.1840 - 0.0000i   0.1938 - 0.0000i   0.2061 - 0.0000i   0.2215 - 0.0000i   0.2408 + 0.0000i   0.2648 + 0.0000i   0.2953 + 0.0000i   0.3341 + 0.0000i   0.3845 - 0.0000i   0.4512 + 0.0000i   0.5416 + 0.0000i   0.6680 + 0.0000i   0.8520 - 0.0000i   1.1338 - 0.0000i   1.5960 + 0.0000i   2.4285 + 0.0000i   4.1447 - 0.0000i   8.4239 - 0.0000i  15.5325 - 0.0000i

syms f(theta) k theta       % Analytic Computation
f(theta) = theta^4 + 1

f(theta) = 

F(k) = int(f*exp(-1j*k*theta), theta, -pi, pi)/(2*pi)

F(k) = 

F = simplify(F, 500)

F(k) = 

kv = sym(1E-50:49);
Fv = double(F(kv))

Fv = 1×50

20.4818 -15.4784 8.3696 -4.0902 2.3737 -1.5407 1.0781 -0.7957 0.6110 -0.4837 0.3924 -0.3246 0.2730 -0.2328 0.2008 -0.1750 0.1538 -0.1363 0.1216 -0.1092 0.0985 -0.0894 0.0815 -0.0745 0.0685 -0.0631 0.0583 -0.0541 0.0503 -0.0469

.

[Paul Anhalt]

Thank you very much! Your solution works even when I make n = 51 to get a size 50 output vector. I don't have an intimate knowledge of FFT. Why does truncating the input vector by 1 eliminate the complex component of the output? When I try:

fft(f, n-2) / (n-2)
fft(f, n-3) / (n-3)
...

The outputs still have imaginary components and stray further from the analytical solution.

[Paul Anhalt]

Interestingly enough, if I leave a small 1-sample gap at the end of my domain, I get a similarly good answer without needing to truncate my vector input. I think this solution uses the same concept as yours where somehow we leave room in the domain for the FFT to make a clean periodic function from a non-periodic input function.

Instead of going from exactly , I go from

n = 50;
f = linspace(-pi, pi - 2 * pi / (n - 1), n);    % Establish domain
f = f.^4 + 1;                                   % Calculate f(theta)
f = fft(f) / n;                                 % FFT
f = f .* (-ones(1, n)).^(0:n - 1);              % Sign pattern of analytical sol'n      

Even though I still get a sizeable imaginary component with this version, my preconditioned conjugate gradient code runs fine. With my first iteration of finding vector f in the original post, the PCG never converged.

Thanks again.

[Star Strider]

As always, my pleasure!

This is an ‘unusual’ Fourier transform, at least in my experience. That likely has more to do with the function being transformed than anything specific to the Fourier transform, or the fft algorithm, however I am not certain for the reason for that result. One reason I did the symbolic transformation was to gain insight into that property, however I did not find that particularly enlightening.

I suspect that the choice of the function to generate the matrix has much to do with the interesting results. The theory underlying that choice may shed light on the result. It would likely require a more detailed investigation into the properties of that function and its Fourier transform to determine that.

Experimenting —

n = 500;
f = linspace(-pi, pi, n);   % Establish domain
f = f.^4 + 1;               % Calculate f(theta)
f = fft(f,(n-1))/(n-1)      % FFT              NOTE: Need To Divide By 'n'

f = 
  20.4823 + 0.0000i  15.4789 + 0.0000i   8.3701 - 0.0000i   4.0907 + 0.0000i   2.3742 + 0.0000i   1.5413 + 0.0000i   1.0786 + 0.0000i   0.7962 + 0.0000i   0.6115 + 0.0000i   0.4843 + 0.0000i   0.3929 - 0.0000i   0.3252 - 0.0000i   0.2735 + 0.0000i   0.2333 + 0.0000i   0.2013 - 0.0000i   0.1755 + 0.0000i   0.1544 - 0.0000i   0.1368 - 0.0000i   0.1221 + 0.0000i   0.1097 - 0.0000i   0.0991 - 0.0000i   0.0899 - 0.0000i   0.0820 - 0.0000i   0.0751 + 0.0000i   0.0690 + 0.0000i   0.0636 - 0.0000i   0.0589 + 0.0000i   0.0546 + 0.0000i   0.0508 + 0.0000i   0.0474 + 0.0000i   0.0444 + 0.0000i   0.0416 + 0.0000i   0.0391 + 0.0000i   0.0368 - 0.0000i   0.0347 + 0.0000i   0.0327 - 0.0000i   0.0310 - 0.0000i   0.0294 + 0.0000i   0.0279 + 0.0000i   0.0265 - 0.0000i   0.0252 + 0.0000i   0.0240 - 0.0000i   0.0229 + 0.0000i   0.0219 + 0.0000i   0.0209 + 0.0000i   0.0200 - 0.0000i   0.0192 + 0.0000i   0.0184 + 0.0000i   0.0177 - 0.0000i   0.0170 - 0.0000i   0.0163 + 0.0000i   0.0157 - 0.0000i   0.0151 + 0.0000i   0.0146 + 0.0000i   0.0141 - 0.0000i   0.0136 + 0.0000i   0.0131 + 0.0000i   0.0127 - 0.0000i   0.0123 + 0.0000i   0.0119 + 0.0000i   0.0115 - 0.0000i   0.0111 + 0.0000i   0.0108 - 0.0000i   0.0105 - 0.0000i   0.0102 + 0.0000i   0.0099 + 0.0000i   0.0096 + 0.0000i   0.0093 + 0.0000i   0.0091 + 0.0000i   0.0088 + 0.0000i   0.0086 + 0.0000i   0.0084 - 0.0000i   0.0082 - 0.0000i   0.0080 - 0.0000i   0.0078 - 0.0000i   0.0076 - 0.0000i   0.0074 - 0.0000i   0.0072 - 0.0000i   0.0070 - 0.0000i   0.0069 - 0.0000i   0.0067 + 0.0000i   0.0066 + 0.0000i   0.0064 - 0.0000i   0.0063 - 0.0000i   0.0061 - 0.0000i   0.0060 - 0.0000i   0.0059 + 0.0000i   0.0058 + 0.0000i   0.0057 - 0.0000i   0.0055 + 0.0000i   0.0054 - 0.0000i   0.0053 + 0.0000i   0.0052 + 0.0000i   0.0051 + 0.0000i   0.0050 + 0.0000i   0.0049 + 0.0000i   0.0048 - 0.0000i   0.0048 - 0.0000i   0.0047 - 0.0000i   0.0046 + 0.0000i   0.0045 + 0.0000i   0.0044 - 0.0000i   0.0044 - 0.0000i   0.0043 + 0.0000i   0.0042 + 0.0000i   0.0042 + 0.0000i   0.0041 + 0.0000i   0.0040 + 0.0000i   0.0040 + 0.0000i   0.0039 + 0.0000i   0.0038 - 0.0000i   0.0038 - 0.0000i   0.0037 - 0.0000i   0.0037 + 0.0000i   0.0036 + 0.0000i   0.0036 + 0.0000i   0.0035 + 0.0000i   0.0035 - 0.0000i   0.0034 + 0.0000i   0.0034 - 0.0000i   0.0033 - 0.0000i   0.0033 - 0.0000i   0.0032 - 0.0000i   0.0032 + 0.0000i   0.0032 - 0.0000i   0.0031 - 0.0000i   0.0031 + 0.0000i   0.0030 - 0.0000i   0.0030 - 0.0000i   0.0030 + 0.0000i   0.0029 + 0.0000i   0.0029 + 0.0000i   0.0029 - 0.0000i   0.0028 + 0.0000i   0.0028 + 0.0000i   0.0028 - 0.0000i   0.0027 + 0.0000i   0.0027 - 0.0000i   0.0027 + 0.0000i   0.0027 + 0.0000i   0.0026 - 0.0000i   0.0026 + 0.0000i   0.0026 - 0.0000i   0.0025 + 0.0000i   0.0025 - 0.0000i   0.0025 + 0.0000i   0.0025 - 0.0000i   0.0025 + 0.0000i   0.0024 + 0.0000i   0.0024 - 0.0000i   0.0024 - 0.0000i   0.0024 - 0.0000i   0.0023 - 0.0000i   0.0023 + 0.0000i   0.0023 - 0.0000i   0.0023 + 0.0000i   0.0023 - 0.0000i   0.0022 - 0.0000i   0.0022 - 0.0000i   0.0022 + 0.0000i   0.0022 + 0.0000i   0.0022 + 0.0000i   0.0022 - 0.0000i   0.0021 + 0.0000i   0.0021 + 0.0000i   0.0021 + 0.0000i   0.0021 + 0.0000i   0.0021 - 0.0000i   0.0021 + 0.0000i   0.0020 + 0.0000i   0.0020 + 0.0000i   0.0020 + 0.0000i   0.0020 + 0.0000i   0.0020 + 0.0000i   0.0020 + 0.0000i   0.0020 - 0.0000i   0.0020 + 0.0000i   0.0019 + 0.0000i   0.0019 + 0.0000i   0.0019 + 0.0000i   0.0019 - 0.0000i   0.0019 - 0.0000i   0.0019 + 0.0000i   0.0019 + 0.0000i   0.0019 + 0.0000i   0.0019 - 0.0000i   0.0018 + 0.0000i   0.0018 - 0.0000i   0.0018 - 0.0000i   0.0018 + 0.0000i   0.0018 - 0.0000i   0.0018 - 0.0000i   0.0018 - 0.0000i   0.0018 - 0.0000i   0.0018 - 0.0000i   0.0018 - 0.0000i   0.0018 - 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 + 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 - 0.0000i   0.0016 + 0.0000i   0.0016 + 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 - 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 + 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0017 - 0.0000i   0.0018 + 0.0000i   0.0018 + 0.0000i   0.0018 + 0.0000i   0.0018 + 0.0000i   0.0018 + 0.0000i   0.0018 + 0.0000i   0.0018 + 0.0000i   0.0018 - 0.0000i   0.0018 + 0.0000i   0.0018 + 0.0000i   0.0018 - 0.0000i   0.0019 + 0.0000i   0.0019 - 0.0000i   0.0019 - 0.0000i   0.0019 - 0.0000i   0.0019 + 0.0000i   0.0019 + 0.0000i   0.0019 - 0.0000i   0.0019 - 0.0000i   0.0019 - 0.0000i   0.0020 - 0.0000i   0.0020 + 0.0000i   0.0020 - 0.0000i   0.0020 - 0.0000i   0.0020 - 0.0000i   0.0020 - 0.0000i   0.0020 - 0.0000i   0.0020 - 0.0000i   0.0021 + 0.0000i   0.0021 + 0.0000i   0.0021 - 0.0000i   0.0021 - 0.0000i   0.0021 - 0.0000i   0.0021 - 0.0000i   0.0022 + 0.0000i   0.0022 - 0.0000i   0.0022 - 0.0000i   0.0022 - 0.0000i   0.0022 + 0.0000i   0.0022 + 0.0000i   0.0023 + 0.0000i   0.0023 - 0.0000i   0.0023 + 0.0000i   0.0023 - 0.0000i   0.0023 + 0.0000i   0.0024 + 0.0000i   0.0024 + 0.0000i   0.0024 + 0.0000i   0.0024 - 0.0000i   0.0025 - 0.0000i   0.0025 + 0.0000i   0.0025 - 0.0000i   0.0025 + 0.0000i   0.0025 - 0.0000i   0.0026 + 0.0000i   0.0026 - 0.0000i   0.0026 + 0.0000i   0.0027 - 0.0000i   0.0027 - 0.0000i   0.0027 + 0.0000i   0.0027 - 0.0000i   0.0028 + 0.0000i   0.0028 - 0.0000i   0.0028 - 0.0000i   0.0029 + 0.0000i   0.0029 + 0.0000i   0.0029 - 0.0000i   0.0030 - 0.0000i   0.0030 + 0.0000i   0.0030 + 0.0000i   0.0031 - 0.0000i   0.0031 + 0.0000i   0.0032 + 0.0000i   0.0032 - 0.0000i   0.0032 + 0.0000i   0.0033 + 0.0000i   0.0033 + 0.0000i   0.0034 + 0.0000i   0.0034 - 0.0000i   0.0035 + 0.0000i   0.0035 - 0.0000i   0.0036 - 0.0000i   0.0036 - 0.0000i   0.0037 - 0.0000i   0.0037 + 0.0000i   0.0038 + 0.0000i   0.0038 + 0.0000i   0.0039 - 0.0000i   0.0040 - 0.0000i   0.0040 - 0.0000i   0.0041 - 0.0000i   0.0042 - 0.0000i   0.0042 - 0.0000i   0.0043 - 0.0000i   0.0044 + 0.0000i   0.0044 + 0.0000i   0.0045 - 0.0000i   0.0046 - 0.0000i   0.0047 + 0.0000i   0.0048 + 0.0000i   0.0048 + 0.0000i   0.0049 - 0.0000i   0.0050 - 0.0000i   0.0051 - 0.0000i   0.0052 - 0.0000i   0.0053 - 0.0000i   0.0054 + 0.0000i   0.0055 - 0.0000i   0.0057 + 0.0000i   0.0058 - 0.0000i   0.0059 - 0.0000i   0.0060 + 0.0000i   0.0061 + 0.0000i   0.0063 + 0.0000i   0.0064 + 0.0000i   0.0066 - 0.0000i   0.0067 - 0.0000i   0.0069 + 0.0000i   0.0070 + 0.0000i   0.0072 + 0.0000i   0.0074 + 0.0000i   0.0076 + 0.0000i   0.0078 + 0.0000i   0.0080 + 0.0000i   0.0082 + 0.0000i   0.0084 + 0.0000i   0.0086 - 0.0000i   0.0088 - 0.0000i   0.0091 - 0.0000i   0.0093 - 0.0000i   0.0096 - 0.0000i   0.0099 - 0.0000i   0.0102 - 0.0000i   0.0105 + 0.0000i   0.0108 + 0.0000i   0.0111 + 0.0000i   0.0115 + 0.0000i   0.0119 - 0.0000i   0.0123 + 0.0000i   0.0127 + 0.0000i   0.0131 - 0.0000i   0.0136 - 0.0000i   0.0141 + 0.0000i   0.0146 - 0.0000i   0.0151 - 0.0000i   0.0157 + 0.0000i   0.0163 - 0.0000i   0.0170 + 0.0000i   0.0177 + 0.0000i   0.0184 - 0.0000i   0.0192 - 0.0000i   0.0200 + 0.0000i   0.0209 - 0.0000i   0.0219 - 0.0000i   0.0229 - 0.0000i   0.0240 + 0.0000i   0.0252 - 0.0000i   0.0265 + 0.0000i   0.0279 - 0.0000i   0.0294 - 0.0000i   0.0310 + 0.0000i   0.0327 + 0.0000i   0.0347 - 0.0000i   0.0368 + 0.0000i   0.0391 - 0.0000i   0.0416 - 0.0000i   0.0444 - 0.0000i   0.0474 - 0.0000i   0.0508 - 0.0000i   0.0546 - 0.0000i   0.0589 - 0.0000i   0.0636 + 0.0000i   0.0690 - 0.0000i   0.0751 - 0.0000i   0.0820 + 0.0000i   0.0899 + 0.0000i   0.0991 + 0.0000i   0.1097 + 0.0000i   0.1221 - 0.0000i   0.1368 + 0.0000i   0.1544 + 0.0000i   0.1755 - 0.0000i   0.2013 + 0.0000i   0.2333 - 0.0000i   0.2735 - 0.0000i   0.3252 + 0.0000i   0.3929 + 0.0000i   0.4843 - 0.0000i   0.6115 - 0.0000i   0.7962 - 0.0000i   1.0786 - 0.0000i   1.5413 + 0.0000i   2.3742 - 0.0000i   4.0907 + 0.0000i   8.3701 + 0.0000i  15.4789 - 0.0000i

So the same result occurs for different values of ‘n’ and the result also tends to converge on the desired magnitudes of the ‘f’ vector for increasing velues of ‘n’.

.

[Paul]

@Paul Anhalt

this comment has the correct concept, but a minor error in implementation. The sampling period and the number of samples of f(theta) have to be selected such that the periodic extension of the sampled sequence would be samples of the periodic extension of the continous signal f(theta). Which is just another way of saying "we leave room in the domain for the FFT to make a clean periodic function from a non-periodic input function." The reason the code in the linked comment still yields a sizable imaginary component is because the sampling of theta is not quite correct.

Code from linked comment:

n = 50;
f = linspace(-pi, pi - 2 * pi / (n - 1), n);    % Establish domain
f = f.^4 + 1;                                   % Calculate f(theta)
f = fft(f) / n;                                 % FFT
figure
subplot(211);stem(0:49,real(f));
subplot(212);stem(0:49,imag(f));

Corrected code

% change the second argument to linspace
f = linspace(-pi, pi - 2 * pi /n, n);           % Establish domain
f = f.^4 + 1;                                   % Calculate f(theta)
f = fft(f) / n;                                 % FFT
figure
subplot(211);stem(0:49,real(f));
subplot(212);stem(0:49,imag(f));

Also, even though n = 50 in this code, that's not the same "n" as in your problem statement. With n = 50 in the code, only the first 25 values in f can be used as approximations to alpha(k), k = 0:24.

[Paul Anhalt]

That's very cool, Paul. Your explanation of how the sampling needs to be done so that the samples of the extended periods line up with each other makes sense, and is a better way to say what was in my head at the time. This thread has given me a lot of intuitive understanding of the FFT.

[Walter Roberson]

It is common for people to build a periodic signal that ends returning to the starting point. For example a triangle that starts and ends with 0. But fft treats the input as indefinitely repeatable, and if you put two cycles of the data together then in the middle you would have the 0 from the end of the signal and then another zero from the beginning, giving you two zeros in a row when there should only be one. Thus you need your signal to end one sample earlier than you would normally expect.

I typically construct the signal as returning, but then throw away the last sample.

Answer 2

[Paul]

Define the function to evaluate f(theta)

ffunc = @(theta) (theta.^4 + 1); % only valid between -pi and pi, zero otherwise

Define a function to evaluate the periodic extension of f(theta) with period 2*pi.

ffuncext = @(theta) ffunc(mod(theta + pi,2*pi) - pi); % periodic extension

We want to use the DFT (via fft()) to compute the discrete Fourier series coefficients of a sequence of samples of f(theta). Those coefficients approximate alpha_k, which are the continuous Fourier series coefficients of f(theta).

In order for this to go smoothly, the sampling period in theta has to be 2*pi/N, where N is an integer. Then we can take N samples of theta and compute the N-point DFT. However, recall that the upper half of the DFT actually corresponds to frequencies from -pi to 0, which correspond to negative values of k. The problem statement only asks for alpha(k) for k = 0:n-1. So we have to choose N to be 2*(n-1), and then only use the first N/2 values of the DFT. For example, let n - 1 = 49

n = 50;
N = 2*(n-1);

The DFT always assumes the first point in the sequence is at index = 0. So, in general, we either sample the function starting from theta = 0 or we can sample starting from theta = -pi and shift the resulting DFT. For this particular f(theta) it won't matter, but I prefer the former.

thetavals = (0:N-1)/N*2*pi; % N samples of theta over one period of f(theta)

Evaluate the function over one period

fvals = ffuncext(thetavals);

Compute the DFT and normalize

c_k = fft(fvals)/N;  % normalized

Compute the exact result

syms f(theta) k
f(theta) = theta^4 + 1;
alpha(k) = simplify(int(f*exp(-1j*k*theta), theta, -pi, pi)/(2*sym(pi)));
alpha(k) = piecewise(k == 0, limit(alpha(k),k,0),alpha(k));

Compare the approximate solution and the exact solution

kvals = 0:n-1;
alphavals = double(alpha(kvals));

figure
stem(kvals(1:n),real(c_k(1:n)))
hold on;
stem(kvals,real(alphavals),'x')

figure
stem(kvals(1:n),imag(c_k(1:n)))
hold on
stem(kvals,imag(alphavals),'x')

Sampling doesn't come for free. The cost is that the c_k only approximate the alpha_k, and the approximation becomes worse for smaller n (fewer samples) and as k inreases towards the Nyquist frequency.

figure
stem(kvals(1:n),abs(c_k(1:n) - alphavals))