am I correct in generating signals and fft?

10 views (last 30 days)
By looking at some examples and searching the internet, I was able to write the following. I don't understand fft resolution. Can we write the number we want? Or how should I write? and, what does "real and imaginary parts of the FFT of x (n)" mean ? Thanks.
Questions asked to me:
* Consider a unit pulse signal of width 65, centered at the origin. This is an even signal and we wish to look into the calculation of its FFT. Generate a signal x(n) of length 128 where the first 33 samples are ones, the next 63 are zeros andthe last 32 samples are ones. Explain how this signal relates to the pulse, andwhat the effect of the zeros is. What is the length of the FFT?
* Calculate the FFT of x (n). Plot the magnitude, real and imaginary parts of the FFT of x (n), and discuss how these results compare with theoretical ones.
and here is what I could write:
clc
clear all
close all
n1=-32:32;
u=[zeros(1,32) 1 zeros(1,32)];
figure("Name", "Ismail Bilgec", "NumberTitle", "off" );
subplot(611);stem(n1,u);title('unit impulse signal');
n2=0:127;
xn=[ones(1,33) zeros(1,63) ones(1,32)];
subplot(612);stem(n2,xn);title('x(n) signal')
X1=fftshift(abs(fft(u,50)/numel(u)));
X2=fftshift(abs(fft(xn,50)/numel(xn)));
w=linspace(-pi,pi,50);
subplot(613);stem(w/pi,X1,'g');title('fft of unit impulse signal');
subplot(614);stem(w/pi,X2,'r');title('fft of x(n) signal');
subplot(615);stem(w, real(X2)); title('real part of the FFT of x(n)');
subplot(616);stem(w, imag(X2)); title('imaginary part of the FFT of x(n)');
length(X1)
length(X2)
  2 Comments
ISMAIL BILGEC
ISMAIL BILGEC on 26 Mar 2022
The problem is to Calculate the FFT of x (n) , not iFFT.

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 25 Mar 2022
The fft result should be normalised by dividing it by the number of elements in the vector:
X1=fftshift(abs(fft(u,50)/numel(u)));
X2=fftshift(abs(fft(xn,50)/numel(xn)));
Other than that, I do not see any problems with the code.
n1=-32:32;
u=[zeros(1,32) 1 zeros(1,32)];
subplot(411);stem(n1,u);title('impulse signal');
n2=0:127;
xn=[ones(1,33) zeros(1,63) ones(1,32)];
subplot(412);stem(n2,xn);title('x(n) signal')
X1=fftshift(abs(fft(u,50)/numel(u)));
X2=fftshift(abs(fft(xn,50)/numel(xn)));
w=linspace(-pi,pi,50);
subplot(413);stem(w/pi,X1,'g');title('fft of unit impulse signal');
subplot(414);stem(w/pi,X2,'r');title('fft of x(n) signal');
length(X1)
ans = 50
length(X2)
ans = 50
.

More Answers (1)

Paul
Paul on 27 Mar 2022
We know that zero-padding in the time domain corresponds to interpolation in the frequency domain. This problem can illustrate the dual property that zero-padding in the frequency domain corresponds to interpolation in the time domain. This duality shouldn't be a surprise as frequency-time duality shows up quite a bit with the DFT. This problem also illustrates other properties of the DFT.
Let's define a continuous time signal
s = @(t) sinc(t);
and evalutate it over -32 < t < 32
tvals = -32:.01:32;
svals = s(tvals);
Now, sample this signal with a sampling period of T = 1 second for 65 samples symmetric around t = 0 to create a 65-point finite duratin sequence u[n]
T = 1;
N = 65;
nvals = -(N-1)/2:(N-1)/2;
uvals = s(nvals*T);
In the discrete-time domain, uvals is the unit pulse sequence defined in the problem
figure
plot(tvals,svals)
hold on
stem(nvals*T,uvals)
xlabel('Time (sec)')
A very important fact that we'll come back to later is that this unit pulse sequence could have been obtained by sampling a different continous-time signal, for example a triangular pulse.
Let's take the DFT of uvals. Recall that the DFT, as implemented by fft(), assumes the first point in the input sequence corresponds to n = 0. But the first point in uvals corresponds to n = -32. So we can either use the time shifting property of the DFT or we can apply ifftshift to uvals.
uvals = ifftshift(uvals);
U = fft(uvals);
From the DFT property "imaginary part in frequency", the imaginary part of U[n] should be zero.
any(uvals - conj([uvals(1) fliplr(uvals(2:end))]))
ans = logical
0
However, U[n] has a small rounding error in its imaginary part, so get rid of it
max(abs(imag(U)))
ans = 1.7356e-31
U = real(U);
Now plot U[n]
figure;
hold on
xline(N/2,'r')
stem(0:N-1,U)
xlabel('Index')
The DFT of the unit pulse sequence is a sequence of shifted unit pulses. The vertical red line corresponds to the highest frequency in discrete time. Points just to the left of the red line correspond to large positive frequencies and points just to the right of the red line correspond to large negative frequencies.
Now, insert 63 zeros into the middle of U[n] and plot the result
X = [U(1:33) zeros(1,63) U(34:end)];
M = numel(X);
figure
hold on
xline(M/2,'r');
stem(0:M-1,X)
xlabel('Index')
Clearly, X[n] is the same as x(n) in the problem. So we can say that X[n] is simply the DFT of the unit pulse with zero-padding at high frequency.
X[n] is still conceptually in the frequency domain, so it is the DFT of some yet-to-be-defined 128-point finite-duration sequence; let's call it x[n] (this is not x(n) in the problem!). To find x[n] we can simply take the iDFT of X[n]
xvals = ifft(X);
We also know from the "imaginary part in time property" of the DFT that x[n] is real so let's verify that and then ensure that against numerical error
any(X - conj([X(1) fliplr(X(2:end))]))
ans = logical
0
max(abs(imag(xvals)))
ans = 0
xvals = real(xvals);
But the problem statement says to take the DFT of X[n], not the iDFT of X[n]. What is the DFT of X[n]?
By the properties of the DFT, we know that for a general discrete-time sequence z[n] of length K that with dft(z[n]) = Z[n]
z = idft(Z) = 1/K*conj(dft(conj(Z)))
For example:
rng(100);
z = rand(1,6); K = 6; Z = fft(z);
[z; ifft(Z) ; 1/K*conj(fft(conj(Z)))]
ans = 3×6
0.5434 0.2784 0.4245 0.8448 0.0047 0.1216 0.5434 0.2784 0.4245 0.8448 0.0047 0.1216 0.5434 0.2784 0.4245 0.8448 0.0047 0.1216
Applying this property to the problem at hand where X[n] = dft(x[n]) and x[n] = idft(X[n]) we have
x = 1/M * conj(dft(conj(X))) = 1/M * conj(dft((X))) because X is real. Therefore, we have
conj(M*x) = M*x = dft(X) because x is real.
So so if we take the DFT of X[n], we get M*x[n]. Verify
norm(M*xvals - fft(X))
ans = 1.7360e-14
But what is x[n]? Recall that this whole process started by applying ifftshift() to uvals and we finished with xvals. So undo that shift
xvals = fftshift(xvals);
As it turns out, x[n], with appropriate scaling by N/M, are points that approximately (maybe exactly) lie on s(t)
NoverM = N/M;
figure
hold on
plot(tvals,svals)
stem(ceil(-M/2:(M/2-1))*T*NoverM,xvals/NoverM);
xlabel('Time (sec)')
figure
hold on
plot(tvals,svals)
stem(ceil(-M/2:(M/2-1))*T*NoverM,xvals/NoverM);
xlim([-10 10])
xlabel('Time (sec)')
The time domain interpolation u[n] is such a good match to s(t) undoubtedly because u[n] resulted from sampling a sinc() function (sinc functions tend to show up a lot in signal processing). But, as stated at the outset, u[n] could have been samples from some other continuous-time signal, like a triangular pulse, in which case the resulting sequence x[n] would be a good approximation to the triangle over the interval -1 < t < 1.
figure
hold on
plot(tvals,triangularPulse(-1,0,1,tvals))
stem(ceil(-M/2:(M/2-1))*T*NoverM,xvals/NoverM);
xlim([-10 10])
xlabel('Time (sec)')
  1 Comment
ISMAIL BILGEC
ISMAIL BILGEC on 27 Mar 2022
Thank you very much for your interest and patience.

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!