DSP Puzzle: What network will create this waveform?

Question

What network consisting of two unit sample delays, three multipliers (with constant gain coefficients), and three adders can we create to convert a square wave into the waveform as shown below:

The waveform is a 1 KHz square wave sampled at 100 KHz.

This is a “DSP Puzzle”, please preface your answer with spoiler notation by typing ythe following two characters first ">!"

Answer 1

A solution with 1 multiply, 2 adds and 2 delays :-)

The target wave form is simply the square wave itself with the fundamental removed. Since the frequency is known, I simply implemented a local oscillator with the same frequency, amplitude and phase of the fundamental and than subtract it out. The most efficient way to create a local oscillator is a simple recursion $$y[n] = 2\cos(\omega_0)\cdot y[n-1] - y[n-2], \quad x_{out}[n] = x_{in}[n]-y[n]$$ The tricky bit is to seed the oscillator states correctly. Assuming we want to implement$$y[n] = A \cdot \cos(\omega_0\cdot n + \varphi) $$ we need to seed $y[n-2] = y[-2]$ and $y[n-1] = y[-1]$. This isn't the most stable oscillator and the amplitude may drift over time. An alternative would be to use a rotating phasor instead, but that would take 4 multiplies.

Results (after fixing a stupid mistake).

Here is the code:

%% DSP Puzzle Dan Boschen 10/1/2023
fs = 100000; % sample rate
f0 = 1000; % square wave frequency
nx = 10000; % number of samples
% create the sqaure wave. This is easy, since it's an integer ratio
np = fs/f0/2;
x0 = [ones(np,1); -1*ones(np,1)];
x0 = repmat(x0,nx/(length(x0)),1);
% find the angle of the fundamental
fx0 = fft(x0);
ax = angle(fx0(1+f0/fs*nx));
% find the magnitude of the fundamental
amp = abs(fx0(1+f0/fs*nx));
A = amp/nx*2;
% build the sine wave recursively: y[n] = a1*y[n-1]-y[n-2]
z2 = A*cos(om0*(-2)+ax); % state y[n-2]
z1 = A*cos(om0*(-1)+ax); % state y[n-1]
a1 = 2*cos(om0); % a1 coefficient. Pole on the unit circle
x2 = zeros(nx,1); % initialize output
for i = 1:nx
z0 = a1*z1 - z2; % recursive sine wave
x2(i) = x0(i)-z0; % subtract from square wave
z2 = z1; z1 = z0; % update state/delays
end
% plot it
xwin = [0.024 0.038];
clf
subplot(2,1,1);
plot(t,x0,'Linewidth',2); set(gca,'xlim',xwin); grid('on');
title('Input');
xlabel('time in s');
ylabel('Amplitude');
subplot(2,1,2);
plot(t,x2,'Linewidth',2); set(gca,'xlim',xwin); grid('on');
title('Output');
xlabel('time in s');
ylabel('Amplitude');
% calcualte fundamental suppression
fx2 = fft(x2);
amp2 = abs(fx2(1+f0/fs*nx));
fprintf('Fundamental supression = %6.2fdB\n', 20*log10(amp2/amp));

Answer 2

The waveform is as Hilmar properly deduced the result of a square wave with the fundamental frequency removed (so alternatively stated the infinite sum of odd harmonics $n$ starting with the 3rd with each harmonic reduced in amplitude by $1/n$.). The original solution I had in mind was to use a second order notch filter as detailed in DSP.SE Post 31028 and the image below to place a notch at the fundamental frequency to be removed:

I adjust the phase between the suppressed fundamental and remaining tones by shifting the center frequency of the notch slightly, which results in the symmetry demonstrated by the plot presented in the question. This comes at the expense of supressing the harmonic, so doesn't have the impressive 260 dB rejection Hilmar achieved as he reported in the comments under the answer (I have 30 dB). It also requires 3 multipliers, 2 delays and 3 adders which isn't as efficient as what Hilmar offered. What it does provide however is stability for longer term operation.

Python code below:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fft as fft
fs = 100e3
tone = 1e3\ n = np.arange(nsamps)
square = sig.square(2*np.pi*tone/fs *n)
# notch filter at tone
r = .998
offset = 0.00003 # 0.0008 for r=.99, 0.00003 for r = .998
wn = 2 * np.pi * tone / fs+ offset
notched1 = sig.lfilter([1, -2 *np.cos(wn), 1], [1, -2 * r * np.cos(wn), r**2], square)