Bandwidth of a signal using MATLAB

I need to find the bandwidth the signal $x = 50 \operatorname{sinc}(125(t-10))^2$.

(Where $\operatorname{sinc}(x)=\dfrac{\sin (\pi x)}{\pi x}$ [normalized] or $\operatorname{sinc}(x)=\dfrac{\sin x}{x}$ [unnormalized])

How would I go about doing this in MATLAB?

Answer

That actually would be a great interview questions. Let's derive it from scratch:

1. $\text{sinc}(x)$: The Fourier Transform of a $\text{sinc}(x)$ is a rectangle. The normalized $\text{sinc}()$ function results in a bandwidth of 1 Hz (from -0.5Hz to 0.5Hz)



2. $t-10$: The $(t-10)$ term is irrelevant to the bandwidth: time shift is equivalent to multiplying with $e^{-j\omega t_0}$, where $t_0$ is time shift, in the frequency domain. That changes the phase but not the amplitude or bandwidth


3. $125*t$: that is a time compression which results in frequency stretching. Since we compress time by 125 the bandwidth increases by 125, i.e. to 125 Hz (from -62.5 Hz to 62.5 Hz)


4. $(\cdot)^2$: Squaring is multiplication in the time domain which corresponds to convolution in the frequency domain. The convolution of two equal rectangles is a triangle. The bandwidth doubles to 250 Hz (from -125 Hz to 125Hz). In general the bandwidth of the multiplication of two arbitrary signals is the sum of the individual bandwidths. That's a simple consequence of the multiplication laws for sine and cosine: you get only sum and difference frequencies.


5. $50*(\cdot)$: that's a simple scale factor that doesn't change the shape of phase or amplitude.

So overall we would expect the spectrum to be of triangular shape (on a double linear plot) going from -125 Hz to +125 Hz with a total bandwidth of 250 Hz.

Writing the Matlab code for this is not as simple as it looks. The code itself is trivial, however you have to pick a sample rate and a proper time axis. In order to do this you need to assess the general shape of the signals first. For example the time domain signal technically goes from $-\infty$ to $+\infty$ and is centered around $t$ = 10s. t = 0:.01:5, really won't do the trick. So in essence you have to do the theoretical analysis before you can write the code (which tends to be a recurring theme in DSP :-))

% 1000 HZ sample rate, since we now the BW is smaller than 500 Hz
fs = 1000; 
% pick a good time length. Because of the time shift we need significantly more

% than 10s. Technically the sinc() has inifite length so we need to extend
% it far enough to cover the vast part of the energy AND we need to start
% at negative times, not at t = 0;
n = 2^16; % that's about 60s. Should be plenty
t = (-n/2:n/2-1)'/fs;   % time axis in seconds, symmetric around t = 0;
x = 50*sinc(125*(t-10)).^2;
% frequency axis:
f0 = fs/n;  % frequency resoution in Hz
faxis = (-n/2:n/2-1)'*f0;   % frequency axis in Hz seconds, symmetric around f = 0;
fx = fft(x);

plot(faxis,circshift(abs(fx),n/2));
title('Spectrum of 50*sinc(125*(t-10)).^2');
xlabel('Frequency in Hz');