How to plot frequency spectrum of a piecewise function in matlab?

Dear friends,

I want to plot the frequency spectrum of this function:

f(t)=1/2*(1+cos(pi*t)) when -1<t<1

otherwise,f(t)=0

I don't know how to do it

Your help would be highly appreciated!

[Paul]

Check rectangularPulse and fourier

Answer 1

[Paul]

Hi Daniel,

Regarding this comment, the use of fftshift and the multiplication by dt both have to do with using the Discrete Fourier Transform (DFT) of the samples of the signal, computed via fft, to approximate the Continuous Time Fourier Transform (CTFT) of the continuous-time signal.

Define the signal

syms t f real
y(t) = 1/2*(1 + cos(sym(pi)*t))*rectangularPulse(-1,1,t);

Compute its CTFT, independent variable f in Hz.

Y(f) = fourier(y(t),t,2*pi*f)

Y(f) = 

Plot the amplitude spectrum

figure;
fplot(abs(Y(f)),[-2 2])
ylim([0 1.1])

We see that even though Y(f) is not bandlimited, it's basically zero for abs(f) > 2. So our sampling frequency should be greater than 4. Choose Fs to be 4 times the approximate bandlimit and compute the corresponding sampling period (sometimes called dt).

Fs = 4*4; Ts = 1/Fs;

Compute samples of the signal between -1 and 1 where it's non-zero. Plot y(t) and its samples.

tval = -1:Ts:1;
yval = 1/2*(1 + cos(pi*tval));

figure;
fplot(y(t),[-1 1])
hold on
stem(tval,yval)
grid

Compute the DFT of the samples of y(t). Because we are only interested in the amplitude spectrum (or is the phase spectrum also of interest?), we can assume the samples start at t = 0 and take the fft directly.

Ydft = fft(yval);
N = numel(Ydft);

The frequency samples corresponding to the DFT samples are

fval = (0:(N-1))/N*Fs;

Plot the amplitude spectrum of the DFT.

figure
stem(fval,abs(Ydft))
axis padded

We see that samples of Ydft corresponding to fval > Fs/2 look like Y(f) for f < 0. This observation is not a coincidence and can be explained by the the fact that the DFT samples are actually frequency domain samples of the Discrete Time Fourier Transform (DTFT), and the DTFT is periodic, and the "central" period of the DTFT, when properly scaled, approximates the CTFT. Hence the use of fftshift to get to samples of the central period of the DTFT.

Note also that the amplitude of the DFT samples are about a factor of 1/Ts larger than the CTFT (easily seen at f = 0). This follows from the fact that the CTFT is an integral, and the DTFT is basically approximating that integral by a rectangular integration, and the width of each recntangle is Ts. Hence the multiplication by Ts.

Apply both operations and compute the corresponding frequency vector.

Ydft = fftshift(Ydft)*Ts;
fval = ceil(linspace(-N/2,N/2-1,N))/N*Fs;

Compare the result to the CTFT.

figure;
fplot(abs(Y(f)),[-2 2])
ylim([0 1.1])
hold on
stem(fval,abs(Ydft))
xlim([-2 2])

As expected, the fftshifted and scaled DFT samples are approximate samples of the CTFT (within the limits of the Nyquist frequency Fs/2).

We can get more frequency domain samples by zero padding the fft

Ydft = fft(yval,128);
N = numel(Ydft);
Ydft = fftshift(Ydft)*Ts;
fval = ceil(linspace(-N/2,N/2-1,N))/N*Fs;

figure;
fplot(abs(Y(f)),[-2 2])
ylim([0 1.1])
hold on
stem(fval,abs(Ydft))
xlim([-2 2])

[Daniel Niu]

Thank you so much, Paul.

[Daniel Niu]

Dear Paul,

I need more help from you!

in the document of fft, They use Y = fft(X);P2 = abs(Y/L);

that is the Y is divided by signal length, It seems it is not necessary to match the Continuous Time Fourier Transform (CTFT) of the continuous-time signal for example in your above code.

I don't know when we need the fft be divided by L?

thank you!

[Paul]

Hi Daniel,

The question asked about the frequency spectrum of a continous-time signal. By my understanding, that would be the CTFT, both amplitude and phase (I only plotted the amplitude). As shown, one can use the DFT to approximate (very well in this example) specific points on that spectrum. Whether or not it was neccessary to match the CTFT depends on what the question really is.

The frequency spectrum of a discrete-time (or sampled) signal is the Discrete Time Fourier Transform (DTFT). Specific points on DTFT can be computed via the DFT. Either the DTFT or the DFT could be an appropriate frequency spectrum in this case.

I guess it all really depends on what the goal actually is.

The DTFT (or DFT) can be divided by L depending on how you want to interpret and use the result. This Answer has more discussion for the case where the underlying signal is assumed to be periodic that may be helpful.

Answer 2

[Askic V]

If this is a homework, you will have to show what have you done so far.

But here are some useful things. First of all, you need to use piecewise function to define your actual function, and then evaluate it for a input array of time samples.

Then you will have to use fft function.

Everything can be found in the documentation.

Here is the start code:

syms t
y = piecewise(t < -1,0,-1 <= t <= 1,1/2*(1+cos(pi*t)),t > 1,0);
a = -2:0.1:2;
z = subs(y,t,a);
figure
plot(a,z)

The rest can be found here:

https://www.mathworks.com/help/matlab/ref/fft.html

[Daniel Niu]

Thank you Askic,

I found something I can not understand. After read some materials, I write the code below.

But why I need to use Y = fftshift(fft(X,n));

in some code, they use fft(X)*dt

What is the difference?

Your help would be highly appreciated!

% syms t
% f = piecewise(-1<t&&t<1, 1/2*(1+cos(pi*t)), 1<t, 0, t<-1, 0);
clear
syms t
f = piecewise(t < -1,0,-1 < t < 1,1/2*(1+cos(pi*t)),t > 1,0)

f = 

%fplot(f)
Fs = 50;         % Sampling frequency
dt = 1/Fs;           % Sampling period
t = -2:dt:2;     % Time vector
L = length(t);      % Signal length
X=double(subs(f,t));
X=fillmissing(X,'linear','SamplePoints',t);

X = 1×201

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

figure;plot(t,X)
title('Cosine Pulse in Time Domain')
xlabel('Time (t)')
ylabel('f(t)')

         


n = 2^nextpow2(L);
%Convert the Gaussian pulse to the frequency domain.
Y = fftshift(fft(X,n));
%Define the frequency domain and plot the unique frequencies.
f = Fs*(-n/2:n/2-1)/n;

P = abs(Y)/n;         % red

figure(2)
plot(f,P) 
title('Cosine Pulse in Frequency Domain')
xlabel('Frequency (f)')
ylabel('|P(f)|')

Answer 3

[Askic V]

This how I would do it:

%%
clear
close all
% 1/2*(1+cos(pi*t))
syms t
y = piecewise(t < -1,0,-1 <= t <= 1,1/2*(1+cos(pi*t)),t > 1,0);

% Perform FFT and look at the amplitude spectrum
Fs = 50;         % Sampling frequency
dt = 1/Fs;       % Sampling period
t_e = -2:dt:2;   % Time vector
z = double(subs(y,t,t_e));

subplot(211)
plot(t_e,z)
title ('Signal in time domain');
xlabel('Time [s]');
L = 2^nextpow2(numel(z)); % 2^n number of samples to do FFT efficiently
Z = fft(z,L);% calculate FFT
% compute two-sided amplitzde spectrum
A2 = abs(Z/L);
% we just want to show one-sided spectrum
A1 = A2(1:L/2+1);
A1(2:end-1) = 2*A1(2:end-1);
f = Fs*(0:(L/2))/L;

subplot(212)
plot(f,A1)
title ('Amplitude spectrum');
xlabel('Frequency [Hz]');

Now, regarding your questions, you don't have fo use fftshift. It is used only if you want to have two-sided representation of the amplitude spectrum with zero frequency component placed in the middle (so called DC component in the signal). In that case frequencies go from -f to f i.e. from negative to positive.

What we want usually is to have one-sided amplitude spectrum with frequencies from 0 fo some f.

Regarding the other question that includes scaling with factor dt = 1/Fs, I must admit I'm not familiar with it when it comes with amplitude (not power) spectrum. These scalings come into play when calculating power spectrum density, and that is something different.

Answer 4

[Askic V]

Hi @Paul, this is a great explanation. I'm not into signal processing and all I know about Fourier transformation I learned at university and in various web pages on the Internet, where I found conflicting information (or at least my understanding is not sufficient).

I also have a few questions and hope this is the right place to ask. I know the following:

• There are continuous F. transform and Fourier series.
• F. transform of a continuous periodic signal in time domain (CTFT) is actually a discrete set of values in frequency domain i.e. Fourier transform of periodic signal is actually a Fourier series. So for periodic signal, only Fourier series coefficients are calculated since its frequency domain is discrete.
• F. transform of continous, time limited signal in time domain is unlimited continous signal in frequency domain i.e. there exist Fourier transformation in form of integral by definition
• If signal is discrete in time domain, then discrete F. transform (DFT) is defined and its representation in frequency domain is periodic signal. Calculating DFT by definition is slow and inefficient, so fast Fourier transformation algorithm is invented (FFT), but it is only a faster way of calculations which yields to the same result as DFT
• When dealing with real world signals in computer, only discrete time values (samples) of the signal is available and stored in memory (because of limited precision), so only DFT (FFT) can be applied.

So my question is the following:

How to do a proper scaling of amplitude in the frequency domain? If one-sided representation is used (only positive frequencies), values should be the same as amplitudes of sin and cos signal components in the time domain. If double-sided representation s used, these amplitudes are 1/2 of originals, because negative frequencies appear.

In your answer you calculated symbolically continuous FT and inspect the freqency content so you can determine the sample freqency that satisfy Nyquist condition. The amplitude was 1. But then amplitude spectrum of the DFT had two samples of amplitudes 16 and 8.xx.

When i used the code example from the documentation:

Z = fft(z,L);% calculate FFT
% compute two-sided amplitzde spectrum
A2 = abs(Z/L);
% we just want to show one-sided spectrum
A1 = A2(1:L/2+1);
A1(2:end-1) = 2*A1(2:end-1);

I got peak at 0.35, which is probably not correct.

I hope you understand my confusion and'll be able to clarify this.

Thank you very much.

[Paul]

Hi @Askic V,

Not sure why starting a new Answer instead of a comment. Anyway ....

Yes, in continuous-time we have the CTFT and the CFS, with the latter applicable to periodic functions or finite duration functions through periodic extension and windowing.

The CTFT of a period signal is not really a discrete set of values. Like any other CTFT its a function of freuqency over the entire real line. The CTFT is an impulse train of Dirac deltas, each scaled by a corresponding exponential CFS coefficient (and perhaps another scale factor depending on the CFS conventions). So, in that sense the CFS coefficients are needed to devine the CTFT of the periodic signal.

For a finite duration signal, the CTFT integral converges for many signals, maybe even all practical signals of interest. However, I'm not sure that every finite duration signal has a convergent CTFT integral. I think the signal still needs to satisfy the Dirichlet conditions, IIRC.

Starting with an absolutely summable discrete-time signal, the Discrete Time Fourier Transform (DTFT) is a continuous, periodic function of frequency. Like the CTFT, some discrete-time signals that aren't absolutely summable and therefore don't have a convergent DTFT can stil be assigned a DTFT. The Discrete Fourier Transform (DFT) is a set of frequency-domain samples of one period of the DTFT. The DFT is not periodic, but can be made into a perodic signal via periodic extension of needed. I don't see the FFT as being some alternative calculation that happens to yield the same result as the DFT. Rather the FFT is means to compute the DFT, i.e., the FFT is not, in and of itself, a transform. See the Description at fft.

I'm not sure why the DTFT can't be applied to a set of samples stored in computer. Like any other continuous-domain signal it can only be computed at a finite set of points, but you can choose as many points as you want.

Personally, I don't like this one-sided, mulitply by two approach. Have to be careful about what's being multiplied, how the result is being presented, how the result is going to be used, and how it's to be interpreted. This forum has countless threads on the topic.

I can't comment on your result without seeing the entirety of your code.

Also, this Answer may be of interest, at least for periodic signals.

[Askic V]

Thank you very much for your time Paul. BTW, my entire code was given in the previous answer (https://www.mathworks.com/matlabcentral/answers/1914810-how-to-plot-frequency-spectrum-of-a-piecewise-function-in-matlab#answer_1174705).

I'll search on this forum for more information.

[Paul]

The code in your answer has L as the length of the DFT after zero-padding. However, on the doc page for fft all of the examples have L as the length of the signal, i.e., not including any zero padding.