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

# Filtering Square Wave with Butterworth Filter using Fourier Coefficients

13 views (last 30 days)

Show older comments

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');

##### 0 Comments

### Answers (1)

Abhimenyu
on 31 May 2024

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:

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!

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!