-1
$\begingroup$

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)}");

The magnitude and phase plots

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?

$\endgroup$
10
  • 2
    $\begingroup$ You apply a rectangular window to your sine waves. That's convolution with a sinc in the frequency domain. You also add a very strange imaginary part to your real sine waves. No idea why, but at this point it's NOT a sine wave anymore $\endgroup$
    – Hilmar
    Nov 30, 2021 at 21:10
  • $\begingroup$ @Hilmar I didn't created a window but instead, I've multiplied all those harmonics with a unit step function so that I can simulate real life scenarios, i.e. dealing with causal signals. Moreover, the unit step function's Fourier transform is not a sinc function as I haven't created something like a boxcar in the time domain. Additionally, I didn't give a imaginary part for the sines but I did for the zero frequency component in order to get zero phase shifts for the frequency components of 0, 1, 2,..., 10 Hz on the right-hand side of the spectrum. $\endgroup$ Nov 30, 2021 at 21:29
  • 3
    $\begingroup$ @Karakoncolos Hilmar is correct: the effect of your use of the unit step is as if you had multiplied the signal with a rectangular window. The reason is because you haven't really used a unit step (which is infinite in duration) but a truncated unit step (which is only len(t) long). Hilmar's other suggestion that the addition of a 2 * pi * i * heaviside(t) to the signal seems arbitrary and probably not what you want to do. Can you please edit your question and add the reasoning for this? $\endgroup$
    – Peter K.
    Nov 30, 2021 at 21:40
  • $\begingroup$ @PeterK. I've already stated the motivation behind adding that zero-frequency imaginary component in the first comment of mine. But, of course, I've added it by thinking about calculations on paper (Unfortunately, the computer is not able to do so.). I've realised that the unit step function is finite in length which in turn justifies Hilmar's statement about the windowing effect and its consequences. $\endgroup$ Nov 30, 2021 at 21:54
  • $\begingroup$ OK! Reopened. Let's see if that helps with getting you an answer. $\endgroup$
    – Peter K.
    Nov 30, 2021 at 22:36

2 Answers 2

0
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ First of all, I haven't start to derive the frequency transfer function $H(f)$ from when the signal that I've introduced into the source code is being perceived by the algorithm. I've just used "According to the source code..." statement for showing the form of the unit impulse response $h(t)$. Other than that, thank you for the notification about the miscalculation of the scaled imaginary unit step function (I am about to edit the related section in the script.) and other comments on the subject. By the way, I've got the answer for the question. $\endgroup$ Dec 2, 2021 at 21:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.