FFT gives wrong answer for a generating function

7 次查看（过去 30 天）
Paul Anhalt
Paul Anhalt 2022-3-18
评论： Paul 2022-3-20
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?
  2 个评论
Jan
Jan 2022-3-18
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 2022-3-20
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2022-3-18
Ran in:
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
.
  6 个评论
显示 4更早的评论
Paul Anhalt
Paul Anhalt 2022-3-20
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
Walter Roberson 2022-3-20
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.

请先登录，再进行评论。

更多回答（1 个）

Paul
Paul 2022-3-18
Ran in:
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))

类别

在 Help Center 和 File Exchange 中查找有关 Bartlett 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by