I am trying to use the rectangular pulse function, but I am not sure how to use it correctly to solve a Fourier problem.

Ts = 1/100;
t = -4:Ts:4-Ts;

% Now looking to create the rectangular pulse

fx = rectangularPulse(t,2);

After this runs, it says invalid number of arguments. I also tried typing in 'rectpuls' because that's intially what I thought the function was, but that didn't work either.

Answer 1

[Star Strider]

It needs to have the beginning and end defined, and an arbitrary argument for the independent variable (here ‘t’) —

syms omega t 
p = rectangularPulse(-4,4,t)

p = 

P = fourier(p,t,omega)

P = 

P = simplify(P, 500)

P = 

figure
fplot(real(P), [-10 10])
grid

.

[Connor Mondock]

Okay thanks and for plotting the fast Fourier transform would I just use the fft function?

[Star Strider]

As always, my pleasure!

Yes.

The fft function requires a double vector, so defining the square pulse would be somethnig like this —

Ts = 1/100;
t = -6:Ts:6;                                            % Time Vector
s = [zeros(1,200) ones(1,801) zeros(1,200)];            % Signal Vector
FTs = fft(s)/numel(s);
Fv = linspace(-2/Ts, 2/Ts, numel(FTs));                 % Negative Nyquist Frequency To Positive Nyquist Frequency

figure
plot(Fv, fftshift(real(FTs)))                           % Double-Sided Fourier Transform
xlabel('Frequency')
ylabel('Amplitude')
xlim([-10 10])

.

[Connor Mondock]

Okay thanks, now how could I use this to evaluate for the discrete fast fourier transform how does that differ from normal fft? Does this involve trapezoidal integration?

[Star Strider]

As always, my pleasure!

This is the discrete fast Fourier transform.

Discrete fft algorithms involve summations, not numerical integration. Most of them use the Cooley–Tukey FFT algorithm and I believe the FFTW suite (see the References in the fft documentation) does as well.

It would be possible to use trapezoidal integration, however that would be fairly slow, simply because of the number of integrations that would have to be performed.

[Star Strider]

To use trapz to transform the square wave pulse —

Ts = 1/100;
t = -6:Ts:6;                                            % Time Vector
s = [zeros(1,200) ones(1,801) zeros(1,200)];            % Signal Vector
% FTs = fft(s)/numel(s);
w = linspace(-pi*3, pi*3, numel(s));                    % Radian Frequency Vector (Must Be Same Length As 's')
f = exp(-1j*w(:)*t) .* s;                               % Fourier Kernel Matrix
F = trapz(w,f);                                         % Fourier Transform

figure
plot(w, real(F),    w, imag(F))
grid
xlabel('$\frac{\omega}{\pi}$', 'Interpreter','latex')
ylabel('Amplitude')
legend('Real(F)','Imag(F)', 'Location','best')

Experiment to get different results.

—————

EDIT — (15 Apr 2022 at 20:40)

Thanks I did it the way you said to but my t needs to be from -4 to 4 and when I use this it says that the matrix dimensions don't agree why is that?

% clear
Ts = 1/100;
t = -4:Ts:4;                                            
s = ones(1,801);            
% FTs = fft(s)/numel(s);
w = linspace(-pi*3, pi*3, numel(s));                    
f = exp(-1j*w(:)*t) .* s;                               
F = trapz(w,f);  

% Fourier Transform
figure
plot(w, real(F),    w, imag(F))
grid
xlabel('$\frac{\omega}{\pi}$', 'Interpreter','latex')
ylabel('Amplitude')
legend('Real(F)','Imag(F)', 'Location','best')

The reason is that ‘s’ and ‘t’ must have the same number of elements.

Defining:

s = ones(1,801);

enforces that.

Here, ‘s’ is a vector of ones and defines a square pulse, however the content of ‘s’ can otherwise be anything so long as the sizes of ‘s’ and ‘t’ match.

.

[Connor Mondock]

Thank you for all the help!

[Star Strider]

As always, my pleasure!

Answer 2

[Paul]

Hi Connor,

This question brings up some interesting aspects of Fourier analysis, and how transforms of the discrete-time samples relate to the Fourier transform of the underlying continuous-time signal.

First, define the continous-time signal and its continuous time Fourier transform (CTFT)

syms t omega_c
s_c(t) = rectangularPulse(-4,4,t);
S_wc(omega_c) = simplify(fourier(s_c(t),t,omega_c))

S_wc(omega_c) = 

The indpendent variable, omega_c, has units of rad/sec (using default Fourier settings in sympref).

Plot the CTFT over the interval from -10 to 10 rad/sec

figure;
subplot(211);fplot(real(S_wc),[-10 10]); ylabel('Real')
subplot(212);fplot(imag(S_wc),[-10 10]); ylabel('Imag')
xlabel('rad/sec')
legend('CTFT')

Now we want to generate samples of s_c(t) for discrete-time analysis. Define an 801-point, symmetric index vector

n = -400:400;

and the associated time vector with sampling period Ts

Ts = 1/100;
t  = n*Ts;

The first and last points of t fall exactly on the leading and trailing edge of s_c(t) and we have to decide what the samples should be on those edges. s_c(t) is an even signal, so it makes sense to choose equal samples on the edges. Let's choose the edge samples to be 1/2 the amplitude of the pulse.

s_n      = 0*t + 1;
s_n(1)   = 1/2;
s_n(end) = 1/2;

With the samples we can compute the discrete-time Fourier transform (DTFT). First we define the continous-time frequency vector over the frequency interval of interest (note the 10 rad/sec is well below the Nyquist frequency of pi/Ts = 314 rad/sec)

wc = linspace(-10,10,1001); % rad/sec

The discrete-time frequency is then

wn = wc*Ts;

wn has units of rad/sample. We can compute the DTFT directly using its definition. To make it easier, define a simple function that maps our mathematical index, n, to the Matlab array index

sfunc = @(n) s_n(n + 401);
dtft1 = 0*wn;
for nval = n
    dtft1 = dtft1 + sfunc(nval)*exp(-1j*wn*nval);
end

Compare to the CTFT by scaling the DTFT by Ts; scale wn by 1/Ts to convert it to rad/sec

figure
subplot(211);
hold on
fplot(real(S_wc),[-10 10]); ylabel('Real')
plot(wn/Ts,real(Ts*dtft1))
xlim([-10 10])

subplot(212);
hold on
fplot(imag(S_wc),[-10 10]); ylabel('Imag')
plot(wn/Ts,imag(Ts*dtft1))
xlim([-10 10])
xlabel('rad/sec')
legend('CTFT','DTFT1')

The DTFT at DC is

Ts*dtft1(wn == 0)

ans = 8

which is exaclty the value of the CTFT at wc = 0, and is a consequence of the selected sample values at the pulse edges.

The DTFT for s_n can also be computed using the function freqz(). freqz() assumes that the first sample in the input signal corresponds to n = 0, In this problem the first sample corresponds to n = -400, so we apply the shift theorem to the output of freqz()

dtft2 = freqz(s_n,1,wn).*exp(1j*400*wn);

figure
subplot(211);
hold on
fplot(real(S_wc),[-10 10]); ylabel('Real')
plot(wn/Ts,real(Ts*dtft1))
plot(wn/Ts,real(Ts*dtft2))
xlim([-10 10])

subplot(212);
hold on
fplot(imag(S_wc),[-10 10]); ylabel('Imag')
plot(wn/Ts,imag(Ts*dtft1))
plot(wn/Ts,imag(Ts*dtft2))
xlim([-10 10])
xlabel('rad/sec')
legend('CTFT','DTFT1','DTFT2')

Next, we can use the Discrete Fourier Transform (DFT) as implemented by the fft() function. The DFT assumes the underlying signal is periodic with period being the lenght of the signal, which our signal isn't. So we'll need to zero-pad to fake out the DFT, so to speak, by making the period look big. Similarly to freqz(), fft() assumes that the first point in the input corresponds to n = 0. Together these issues can be addressed in a couple of ways. One approach used here is to zero pad the signal on both sides, and then use ifftshift before taking fft()

dft = fft(ifftshift([zeros(1,200) s_n zeros(1,200)]));

The discete-time discrete frequency vector corresponding to the DFT samples is

wdft = (0:numel(dft)-1)/numel(dft)*2*pi; % rad/sample

Note that wdft is one sample short of going all the way around the unit circle, i.e,. the last point is is not 2*pi.

We want our plot to cover the negative frequencies so use fftshift()

dft = fftshift(dft);

and correspondingly modify the frequency vector (there are probably better ways to do this, but I think this way makes it clear what's happening)

wdft = [wdft(wdft > pi-eps(pi)) - 2*pi, wdft(wdft <= pi-eps(pi))];

Because numel(dft) is odd, the first and last points of wdft are

wdft([1 end])

ans = 1×2

-3.1390 3.1390

which are symmetric around 0 and less than pi in absolute value. If nume(dft) were even, then wdft(1) == pi, and wdft(end) would be pi - dwdft, where dwdft is the spaicing between samples of wdft.

The dft should be samples of the DTFT

figure
subplot(211);
hold on
fplot(real(S_wc),[-10 10]); ylabel('Real')
plot(wn/Ts,real(Ts*dtft1))
plot(wn/Ts,real(Ts*dtft2))
stem(wdft/Ts,real(Ts*dft))
xlim([-10 10])

subplot(212);
hold on
fplot(imag(S_wc),[-10 10]); ylabel('Imag')
plot(wn/Ts,imag(Ts*dtft1))
plot(wn/Ts,imag(Ts*dtft2))
stem(wdft/Ts,imag(Ts*dft))
xlim([-10 10])
xlabel('rad/sec')
legend('CTFT','DTFT1','DTFT2','DFT')

Finally, the question asked about the function rectpuls(). Generate samples using rectpuls()

s_nrect = rectpuls(t,8);  

The samples on the edges are

s_nrect([1 end])

ans = 1×2

1 0

So rectpuls() takes the full amplitude at the leading edge and zero at the trailing edge. s_nrect is not symmetric around the origin. Consequently, its DTFT (and DFT) will pick up a small imaginary part and distort the real part relative to the CTFT, though the value at DC is maintained

dtftrect = freqz(s_nrect,1,wn).*exp(1j*400*wn);
Ts*dtftrect(wn == 0)

ans = 8

figure
subplot(211);
hold on
fplot(real(S_wc),[-10 10]); ylabel('Real')
plot(wn/Ts,real(Ts*dtft1))
plot(wn/Ts,real(Ts*dtft2))
stem(wdft/Ts,real(Ts*dft))
plot(wn/Ts,real(Ts*dtftrect))
xlim([-10 10])

subplot(212);
hold on
fplot(imag(S_wc),[-10 10]); ylabel('Imag')
plot(wn/Ts,imag(Ts*dtft1))
plot(wn/Ts,imag(Ts*dtft2))
stem(wdft/Ts,imag(Ts*dft))
plot(wn/Ts,imag(Ts*dtftrect))
xlim([-10 10])
xlabel('rad/sec')
legend('CTFT','DTFT1','DTFT2','DFT','DTFTRect')

[Connor Mondock]

Thanks for all the help, but could someone show me how I could solve a problem like this by numerical integration and/or using the trapezoidal rule?

[Paul]

What exactly do you want to solve via numerical integration?

[Connor Mondock]

creating a vector with a given frequency and then solving it using numerical integration or more specifically using the trapezoidal rule.

[Star Strider]

@Connor Mondock — See my most recent Comment to my answer.

[Connor Mondock]

Thanks I did it the way you said to but my t needs to be from -4 to 4 and when I use this it says that the matrix dimensions don't agree why is that?

clear
Ts = 1/100;
t = -4:Ts:4;                                            
s = [zeros(1,200) ones(1,801) zeros(1,200)];            
% FTs = fft(s)/numel(s);
w = linspace(-pi*3, pi*3, numel(s));                    
f = exp(-1j*w(:)*t) .* s;                               
F = trapz(w,f);  

% Fourier Transform
figure
plot(w, real(F),    w, imag(F))
grid
xlabel('$\frac{\omega}{\pi}$', 'Interpreter','latex')
ylabel('Amplitude')
legend('Real(F)','Imag(F)', 'Location','best')

[Paul]

I think you meant to post this comment in StarStrider's answer.

I can see a few problems, assuming that you're trying to approximate the CTFT.

s and t must both be row vectors with the same number of elements.

The independent variable of integration in the CTFT is t, not w.

The frequency variable, w, has units of rad/sec. Because the plot is against w, the axis label should just be w.

clear
Ts = 1/100;
t = -4:Ts:4;                                            
s = ones(1,801); % just cover the nonzero part of the pulse           
% FTs = fft(s)/numel(s);
w = linspace(-pi*3, pi*3, numel(s));                    
f = exp(-1j*w(:)*t) .* s;                               
F = trapz(t,f,2); % ingerate wrt t, result is a function of w  
% Fourier Transform
figure
plot(w, real(F),    w, imag(F))
grid
xlabel('\omega (rad/sec)')
ylabel('Amplitude')
legend('Real(F)','Imag(F)', 'Location','best')

Result is an essentially identical plot to the others in this answer.