Filtering Square Wave with Butterworth Filter using Fourier Coefficients

Hello there, I'm trying to filter a square wave with a 2nd order butterworth HPF in order to retrieve the fundamental sine wave with minimal attenuation. When I run my code to display the magnitude spectrum of the input and output coefficients, there is an attenuation of about 1.7dB. Yet when I reconstruct the signal, the resulting sine wave appears to be attenuated much more (about half the size of the input, about 6dB). This leads me to believe there is an issue with the second for loop but I am unable to find it. Any help would be much appreciated.

clear; clc;

T = 0.1; %in ms
F = 1/T; %in kHz
Fs = 1000;
offset = 0;
amp = 1;
duty = 50; %percent duty cycle
t = 0:0.001:0.5;
fc = 12; % in kHz

sq_wave = offset+amp*square(2*pi*F.*t,duty); %general square wave function

N=10; %no. of harmonics
nvec = -N:N;
c_in = zeros(size(nvec));
f = nvec*F;

for n = nvec
    m = n+N+1;
    c_in(m) = (1/2)*((sin(pi*n/2))/(pi*n/2)); %sinc function
    
    if (n == 0)
      c_in(m) = 0.0; %DC Offset
    end
end

w = f/fc;
Hf = 1 ./ ((w*1i).^2 + 1.414*(w*1i) + 1);
c_out = c_in .* Hf;

stem(f,abs(c_in),'r','LineWidth',1.5);
hold on
stem(f,abs(c_out),'b','LineWidth',1.5);
hold off
title('Magnitude Spectrum of Filter Output and Input')

A = zeros(2*N+1,length(t));
for n = nvec
    m=n+N+1;
    A(m,:) = c_out(m) .* exp(1i*2*pi*n*F.*t);
end
sine_wave = sum(A);
%plot(t,real(sine_wave),'b',t,sq_wave,'r');

Answer 1

[Abhimenyu]

Hi,

I understand that you have observed difference in the attenuation of the reconstructed sine wave. When you observe a certain attenuation in the magnitude spectrum, it directly relates to the amplitude of the frequency components. However, the perceived discrepancy in the time-domain reconstruction (i.e., the sine wave appearing to be attenuated more than expected) might be due to how these components add up in time, especially considering phase shifts introduced by the filter. Please follow the below-mentioned reasons and recommendations for improving the results:

• Harmonic Content: Given that a square wave is comprised solely of odd harmonics of the fundamental frequency, ensure your analysis meticulously calculates each harmonic's contribution.
• Filter Phase Response: The current reconstruction overlooks the filter's phase response. This aspect, while not influencing the magnitude spectrum, significantly impacts the time-domain signal. The phase shifts could be the key to understanding the perceived extra attenuation.
• Numerical Precision and Aliasing: Confirm that the sampling rate ((F_s)) is sufficiently high to capture the harmonics and filter effects accurately.

Please follow the below-mentioned example MATLAB R2024a code that uses the above recommendations to successfully reconstruct the sine wave:

clear; clc;

% Simulation parameters
T = 0.1; % Period in ms (10 kHz fundamental frequency)
F = 1/T; % Fundamental frequency in kHz
Fs = 1000; % Sampling frequency (kHz)
t = 0:1/Fs:0.5; % Time vector (s)

Cutoff Frequency Adjustment: The cutoff frequency (fc) has been slightly increased to ensure it is above the fundamental frequency. This will more effectively attenuate the fundamental frequency component and allow for a better examination of the high-pass filter's effect.

fc = 15; % Adjusted cutoff frequency in kHz for the HPF, above the fundamental

% Square wave generation
sq_wave = square(2*pi*F.*t, 50); % 50% duty cycle square wave

Odd Harmonics Only: This code now correctly focuses on the odd harmonics of the square wave, which are the only ones present due to its symmetry and Fourier series representation.

Fourier Coefficients for Square Wave: The calculation of Fourier series coefficients (c_in) now directly uses the formula for the amplitudes of the odd harmonics of a square wave, which is (4/pi) * (1/n) for the nth harmonic.

% Fourier series coefficients for square wave, focusing on odd harmonics
nvec = 1:2:(2*10+1); % Considering only odd harmonics up to the 21st
c_in = zeros(size(nvec));
f = nvec*F; % Frequency vector for harmonics

for k = 1:length(nvec)
    n = nvec(k);
    c_in(k) = (4/pi) * (1/n); % Amplitude of nth harmonic
end

% Filter design (2nd order Butterworth HPF)
w = f/fc;
Hf = 1 ./ ((w*1i).^2 + 1.414*(w*1i) + 1); % HPF transfer function

% Apply filter to the coefficients
c_out = c_in .* Hf;

Magnitude Spectrum and Reconstructed Signal: The plotting sections clearly differentiate between the input and output coefficients' magnitude spectra, and the reconstructed signal versus the original square wave.

% Plot magnitude spectrum of filter output and input
figure;
stem(f,abs(c_in),'r','LineWidth',1.5);
hold on;
stem(f,abs(c_out),'b','LineWidth',1.5);
hold off;
title('Magnitude Spectrum of Filter Output and Input');
xlabel('Frequency (kHz)');
ylabel('Magnitude');
legend('Input Coefficients','Output Coefficients');

% Reconstruct the filtered signal
sine_wave = zeros(size(t));
for k = 1:length(nvec)
    n = nvec(k);
    sine_wave = sine_wave + real(c_out(k) * exp(1i*2*pi*n*F.*t));
end

Normalization: After reconstructing the signal, it's normalized to match the peak amplitude of the original square wave. This step ensures that the reconstructed signal does not artificially exceed the original signal's amplitude due to the summation of harmonic components.

% Normalize the reconstructed signal to match the original square wave amplitude
sine_wave = sine_wave / max(abs(sine_wave)) * max(abs(sq_wave));

% Plot the reconstructed signal and original square wave
figure;
plot(t, sine_wave, 'b', 'LineWidth', 1.5);
hold on;
plot(t, sq_wave, 'r--', 'LineWidth', 1.5);
hold off;
title('Reconstructed Signal and Original Square Wave');
xlabel('Time (ms)');
ylabel('Amplitude');
legend('Reconstructed Sine Wave', 'Original Square Wave');

I hope this helps!