How to plot frequency spectrum of a piecewise function in matlab?
23 views (last 30 days)
Show older comments
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!
Accepted Answer
Paul
on 18 Feb 2023
Edited: Paul
on 19 Feb 2023
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)
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])
3 Comments
Paul
on 21 Feb 2023
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.
More Answers (3)
Askic V
on 18 Feb 2023
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
Askic V
on 18 Feb 2023
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.
0 Comments
Askic V
on 20 Feb 2023
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.
3 Comments
Askic V
on 21 Feb 2023
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.
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!