How to plot DSP signals with index values
5 views (last 30 days)
Show older comments
I wrote a code for being able to plot dtft values' magnitude spectrums against their dft index k values. I'm supposed to recreate these graphs, but I can't even get the first graph to look correctly. I managed to get this to work for f/fs as the x axis and I thought all I would have to do is alter the x-axis to make it for k values (N/2). I'm not sure why the curve itself is not plotting correctly or why the signals (vertical lines) are going up to 32 when theoretically they should be matching where the curve ends. Any help or insight on where altering my code went wrong, would be much appreciated.
Edit to add code:
f1 = 2e3;
f2 = 2.5e3;
f3 = 3e3;
fs = 10e3;
t = 0:1/fs:0.1;
x = cos(2*pi*f1*t) + cos(2*pi*f2*t) + cos(2*pi*f3*t);
L = 100;
N1 = 32;
N2 = 64;
w = rectwin(L)';
X1 = fft(x(1:L).*w, 256);
X2 = fft(x(1:L).*w, 256);
X1_dtft = abs(fft(x(1:L).*w, N1));
X2_dtft = fft(x(1:L).*w, N2);
figure(1);
subplot(2,2,1);
plot(0:255, abs(X1), '--');
hold on
stem(0:N1/2-1, abs(X1_dtft(1:N1/2)), 'r', 'linewidth', 0.5, 'marker', 'none')
title('Rectangular window, L=100, N=32');
xlabel('DFT index k');
ylabel('Magnitude Spectrum');
xlim([0 16]);
ylim([0 50]);
3 Comments
Paul
on 4 Apr 2023
Hi Kristina,
Looking just at the code without knowing exactly the problem that the code is trying solve .....
We start with a signal that is sampled at 10 kHz for a fairly long time interval.
Then, we window that signal down to a finite duration signal of 100 points. I'm going to call that signal x. It looks like, in general, the goal is compare the DTFT of x with the DFT of x or possibly frequency samples of the DTFT.
One way to approximate the DTFT of x is to use fft but with a lot of zero padding. That looks like what you're doing with X1 using 256 points. So X1 could be viewed as the DTFT of x. But then X2 is computed the same way; what is supposed to be the difference between X1 and X2? Then, two other variables called X1_dtft and X2_dtft are computed, but those computations aren't the DTFT of x; maybe they are misnamed.
X1_dtft and X2_dtft are computed using fft lengths of N1 = 32 and N2 = 64, which are both smaller than the length of x, which is L = 100. Calling fft with N1 < L truncates the input to the first N1 points and takes the fft of that truncated input. Same with N2. I wonder if that's really the goal. If the goal is get 32- and 64-point samples of the DTFT of x, then it would be more appropriate to use freqz. But it's hard to say w/o more clarification on what the goal is with N1 and N2.
Fundamentally, the independent variable of the DTFT is frequency, whereas for the DFT it is index. For comparing the two, I think it's more typical to convert the DFT indices to frequency samples, though I suppose it's feasible to go the other way and plot the DTFT against a "mapping" from frequency to index.
Accepted Answer
Joseph
on 4 Apr 2023
I think you are running into a problem with frequency alignment between your two signals. The indices corresponding to your DTFT plot are not the same size in frequency as DFT indices. It may be easiest to convert the DTFT to frequency and then scale it back by the DFT step size.
This code is close to working but needs couple more edits
clc
clear all
close all
f1 = 2e3;
f2 = 2.5e3;
f3 = 3e3;
fs = 10e3;
t = 0:1/fs:0.1;
x = cos(2*pi*f1*t) + cos(2*pi*f2*t) + cos(2*pi*f3*t);
L = 100;
N1 = 32;
N2 = 64;
% w = rectwin(L)';
w = ones(1,L);
X1 = fft(x(1:L).*w, 256);
X2 = fft(x(1:L).*w, 256);
X1_dtft = abs(fft(x(1:L).*w, N1));
X2_dtft = fft(x(1:L).*w, N2);
figure(1);
% subplot(2,2,1);
f1 = fftfreq(t(1:L),256); %frequencies
f2 = fftfreq(t(1:L),N1);
DFT_index_step = f2(2)-f2(1);
x1 = f1/DFT_index_step-min(f1/DFT_index_step);
x2 = f2/DFT_index_step-min(f2/DFT_index_step);
plot(x1, abs(X1), '--');
hold on
stem(x2(1:N1), abs(X1_dtft(1:N1)), 'r', 'linewidth', 0.5, 'marker', 'none')
title('Rectangular window, L=100, N=32');
xlabel('DFT index k');
ylabel('Magnitude Spectrum');
xlim([0 16]);
ylim([0 50]);
function f = fftfreq(t,N)
if ~exist('N','var') || isempty(N)
N = length(t);
end
Fs = diff(t);
Fs = 1/Fs(1);
df = Fs/N;
f = -Fs/2:df:(Fs/2-df);
end
0 Comments
More Answers (0)
See Also
Categories
Find more on Fourier Analysis and Filtering 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!