Unexpected Sidebands on the Frequency Spectrum

Question

I am working on an alternative approach for a specific signal processing sequence which involves a kind of a test that I am using for justify the mathematical background of the design. The test is comprised of spreading ten harmonics of a sine function of frequency of 1 Hz over the frequency spectrum and interpreting the resultant magnitude and phase plots. While doing so, in contrast to the theoretical analysis, I’ve observed sinc-like sidebands and, as a result, a very confusing phase graph. The relevant MATLAB source code and consequent graphs are given below:

t = -10: 0.001: 10; 

Ts = 0.001;  % Sampling period in seconds
fs = 1000;  % Sampling frequency in Hertz

function_0 = sin(2 * pi * t) + sin(4 * pi * t) + sin(6 * pi * t) + sin(8 * pi * t) + sin(10 * 
pi * t) + sin(12 * pi * t) + sin(14 * pi * t) + sin(16 * pi * t) + sin(18 * pi * t) + sin(20 * 
pi * t);

function_1 = (2 * pi * i * heaviside(t)) + (4 * pi * function_0 .* heaviside(t));

number_of_samples_1 = length(function_1);
number_of_samples_2 = power(2, nextpow2(number_of_samples_1));

f = linspace(-fs / 2, fs / 2, number_of_samples_2);   

function_2 = fft(function_1, number_of_samples_2) * Ts;
function_3 = fftshift(function_2);
function_4 = abs(function_3);

subplot(2, 1, 1);

plot(f, function_4, "r");
xlabel("Frequency (Hz)");
ylabel("|Function-4 (f)|");

subplot(2, 1, 2);

theta = (180 * angle(function_3) / pi);

plot(f, theta, "r");
xlabel("Frequency (Hz)");
ylabel("\angle{Function-4(f)}");

According to the source code, the unit impulse response is;

$$h(t) = {j2\pi u(t)} + {4\pi \sum_{k=1}^{10} sin(2\pi kt)u(t)}$$

The reason for multiplying the imaginary part of $h(t)$ is achieving the same magnitude in the frequency domain.

We can split the unit impulse response $h(t)$ into two functions as in the expression below:

$$h(t) = a(t) + b(t)$$

The act of taking Fourier transform of $h(t)$ can be modelled as;

$$h(t) \xrightarrow{\mathscr{F}} H(f)$$

By assuming that the system of concern is LTI;

$$a(t) + b(t) \xrightarrow{\mathscr{F}} A(f) + B(f)$$

Let's find Fourier transforms of each function.

$$a(t) = j2\pi u(t) \xrightarrow{\mathscr{F}} A(f) = \frac{j2\pi}{j\omega} + {j2\pi \frac{\delta(f)}{2}}= \frac{1}{f} + j\pi \delta(f)$$

$$b(t) = {4\pi \sum_{k=1}^{10} sin(2\pi kt)u(t)} \xrightarrow{\mathscr{F}} B(f) = 4\pi [\frac{1}{j\omega} \ast \sum_{k=1}^{10} \frac{\delta(f+k) - \delta(f-k)}{j2}] = {\sum_{k=1}^{10} \frac{1}{f-k}} - {\sum_{k=1}^{10} \frac{1}{f+k}}$$

As a result, the frequency transfer function $H(f)$ becomes;

$$H(f) = {\frac{1}{f}} + {j\pi \delta(f)} + {\sum_{k=1}^{10} \frac{1}{f-k}} - {\sum_{k=1}^{10} \frac{1}{f+k}}$$

By creating the causal imaginary DC signal $a(t)$ and using causal sines, the following statements are true for the center and righ-hand side frequency components which are the only components of $H(f)$ that can be observed in real life spectrum measurements (which points out the single-sided Fourier transform);

1. Each component has nearly the same magnitude.

2. Each component's phase shift will be zero which corresponds to zero time delay in the time domain.

Note that the analysis is built upon CTFT and by the assumption that the time interval is not bounded.

What is the cause of those phenomena?

Answer 1

To remove some of the effects of the rectangular window created by using a finite length FFT, try instead doing a DFT whose length is an exact integer multiple of all the input sinusoids. You may have to adjust Ts and your number of samples to make that possible.

Answer 2

There are few things happening here

According to the source code, the unit impulse response is...

Not really. Once you have source code you are operating in discrete time while all your math is in continuous time. In order to make the transition, you need to sample at some finite sample rate. Turns out you cannot sample a causal signal without aliasing. Any causal signal is infinite in frequency so it's not bandlimited.

As soon as you apply an FFT you are also sampling in frequency as well (with the sampling interval determined by the FFT length). Signals that are discrete in one domain are periodic in the other domain, so the whole concept of causality becomes iffy.

$a(t) = j2\pi u(t) \xrightarrow{\mathscr{F}} A(f) = \frac{j2\pi}{j\omega} = \frac{1}{f}$

Nope

$$ \mathcal{F}(u(t)) = \frac{1}{j\omega} + \pi\delta(\omega)$$

You are missing the delta at $\omega = 0$

Overall the concept of a unit step doesn't work well in the context of the discrete Fourier transform (DFT). Since both time and frequency domain signals are periodic a definition of $u[n] = 0, n < 0$ doesn't really work. You could try to zero out the second half of the time domain period but that just gives you a rectangular window again.

function_1 = (2 * pi * i * heaviside(t)) + (4 * pi * function_0 .* heaviside(t));

You are trying to make things causal, but you don't. Since the time domain signal is discrete, zeroing out the "negative" part and then applying FFT shift is exactly equivalent to not time shifting and zeroing out the second half of the time domain buffer. It's a rectangular window again.

Once you apply a unit step to a sine wave, it's not a sine wave anymore, so of course this will give some spectral widening.

Any discrete unit step is problematic. It never decays and if you want to apply an FFT you will have to truncate it. You can implement it as a simple integrator i.e. $y[n] = x[n] + y[n-1]$ but this only marginally stable and the z transform does not converge.