绘制符号傅里叶级数：现有代码疑问及后续开发技术咨询

Hey there! Let's work through fixing your Fourier Series code and getting that plot up and running. First, I spotted a couple of critical bugs in your current implementation that we need to address before diving into plotting:

Key Issues in Your Current Code #

1. Incorrect Coefficient Storage
   In your first for loop, you're adding each computed coefficient to the entire cosCoefficients/sinCoefficients array instead of assigning it to the specific index n. This means all elements in the array will end up with the same (summed) value, which is not what we want.

2. Fourier Series Construction Error
   Your second for loop resets FS to 0.5*a0 + ... on every iteration, so only the final term of the series will be preserved. We need to accumulate each term instead of overwriting the result.

Fixed Fourier Series Function #

Here's the corrected version of your function, with comments explaining each change:

function FS = FourierSeries(f, degree)
    cosCoefficients = zeros(1, degree);
    sinCoefficients = zeros(1, degree);
    syms x;
    
    % Calculate the DC component (a0)
    a0 = double((1/pi)*int(f, -pi, pi));
    
    % Compute cosine and sine coefficients for each harmonic n
    for n = 1:degree
        % Assign to the nth index instead of adding to the whole array
        cosCoefficients(n) = double((1/pi)*int(f*cos(n*x), -pi, pi));
        sinCoefficients(n) = double((1/pi)*int(f*sin(n*x), -pi, pi));
    end
    
    % Build the Fourier Series by accumulating each term
    FS = 0.5*a0; % Start with the DC component
    for n = 1:degree
        FS = FS + cosCoefficients(n)*cos(n*x) + sinCoefficients(n)*sin(n*x);
    end
end

Plotting the Symbolic Fourier Series #

To plot the symbolic expression returned by FourierSeries, we'll convert it to a numerical function handle (using matlabFunction) so we can compute values across a range of x. Here's a complete example using a square wave as the test function:

% Define a test function (periodic square wave, period 2π)
syms x;
f = piecewise(x < 0, -1, x >= 0, 1);

% Compute Fourier Series up to degree 10
degree = 10;
FS = FourierSeries(f, degree);

% Convert symbolic expression to a numerical function for plotting
fs_numeric = matlabFunction(FS);
x_vals = linspace(-2*pi, 2*pi, 1000); % Range to plot over
fs_vals = fs_numeric(x_vals);
original_vals = double(subs(f, x, x_vals)); % Compute original function values

% Create the plot
figure;
plot(x_vals, original_vals, 'k--', 'LineWidth', 1.5, 'DisplayName', 'Original Square Wave');
hold on;
plot(x_vals, fs_vals, 'b-', 'LineWidth', 1, 'DisplayName', ['Fourier Series (Degree ', num2str(degree), ')']);

% Add plot labels and formatting
xlabel('x');
ylabel('f(x)');
title('Fourier Series Approximation of a Square Wave');
legend;
grid on;

Additional Notes #

• If your target function isn't periodic with period 2π, you'll need to adjust the integral bounds and coefficient formulas to match its period.
• Increasing degree will improve the approximation accuracy but will also increase computation time (since we're computing more integrals).
• If MATLAB can't compute an analytical integral for your function, it will fall back to numerical integration automatically, and the double() conversion will still work.

内容的提问来源于stack exchange，提问作者Skm