Spectra replica after sampling a continuous time signal

Please I need to generate and plot the spectra replicas after sampling a continuous time signal , I generate a signal ( time limited to 0.3 seconds ) , composed of 3 different frequencies using samplling frequency fs = 40000 , then I downsample this signal by 10 ( ratio of fs/Fs) , then computed the DFT , but the plot for spectra replicas is not clear at multiples of 4000 Hz , also i didn't get the original spectrum whis is centered at zero , I used hold on to get it but I think this is incoorect , please advise , thanks .

%% Spectra replica after sampling a continuous time signal 
clear
close all
clc
fs=40000; % sampling frequency to generate the original continuous time signal
Fs=4000; % sampling frequency to sample the original signal
dt=1/fs;
Ts=1/Fs;
t=0:dt:0.2-dt; % time vector , duration of the signal = 0.2 seconds 
sig=2*cos(2*pi*400*t)+cos(2*pi*800*t)-3*sin(2*pi*1200*t); % original continuous time signal
Nfft=2^nextpow2(length(sig)); % fft points
Sf=fft(sig,Nfft); % DFT using fft
Sf = fftshift(Sf)/length(sig);   % Shift zero-frequency component to center of spectrum
f= (-Nfft/2:Nfft/2-1)*(fs/Nfft); % zero-centered frequency range
sig_sampled = downsample(sig,fs/Fs);
Sf1=fft(sig_sampled,Nfft); % DFT using fft
Sf1 = fftshift(Sf1)/length(sig_sampled);   % Shift zero-frequency component to center of spectrum
figure(1)
subplot(211)
plot(f/1000,abs(Sf),'LineWidth',1.25) ; grid ; title('Double-Sided Magnitude Spectrum of the baseband signal m(t)') ; 
xlabel('f (KHz)') ; ylabel('|M(f)|') ;  ylim([0 1.6]) ;
subplot(212)
plot(f/1000,abs(Sf),'b','LineWidth',1.25) ; grid ; title('Double-Sided Magnitude Spectrum of the baseband signal m(t)') ; hold on
plot(f/1000,abs(Sf1),'b','LineWidth',1.25) ; grid ; title('Double-Sided Magnitude Spectrum of the sampled signal m(t)') ; ylim([0 1.8])
xlabel('f (KHz)') ; ylabel('|M1(f)|'); 

Answer 1

[Bjorn Gustavsson]

But you do get the original spectra. When you downsample to an effective sampling-frequency of 4000 Hz, your Nyquist-frequency goes down to 2000 Hz, and you'll have to settle for that frequency-range. The peaks separated by multiples of 4000 Hz might have something to do with the discontinuity between the end of the down-sampled signal and the beginning - the FFT is after all a circular Fourier-transform.

HTH

[Mohamad]

I didn't get the original spectrum in the spectrum of the sampled signal S1(f) , I made it in the plot using hold on .

So I don't understand why I didn't get the original spectrum S(f) ?

[Bjorn Gustavsson]

Ah, I see. You'll have to adjust the frequency-vector to account for your reduced sampling-frequency. Try this:

%% Spectra replica after sampling a continuous time signal 
fs = 40000; % sampling frequency to generate the original continuous time signal
Fs = 4000; % sampling frequency to sample the original signal
dt = 1/fs;
Ts = 1/Fs;
t = 0:dt:0.2-dt; % time vector , duration of the signal  =  0.2 seconds 
sig = 2*cos(2*pi*400*t)+cos(2*pi*800*t)-3*sin(2*pi*1200*t); % original continuous time signal
Nfft = 2^nextpow2(length(sig)); % fft points
Sf = fft(sig,Nfft); % DFT using fft
Sf  =  fftshift(Sf)/length(sig);   % Shift zero-frequency component to center of spectrum
f =  (-Nfft/2:Nfft/2-1)*(fs/Nfft); % zero-centered frequency range
sig_sampled  =  downsample(sig,fs/Fs);
Nfft_ds = 2^nextpow2(length(sig));
f_ds =  (-Nfft/2:Nfft/2-1)*((Fs)/Nfft); % zero-centered frequency range
Sf1 = fft(sig_sampled,Nfft_ds); % DFT using fft
Sf1 = fftshift(Sf1)/length(sig_sampled);   % Shift zero-frequency component to center of spectrum
figure(1)
subplot(211)
plot(f/1000,abs(Sf),'LineWidth',1.25) ; grid ; title('Double-Sided Magnitude Spectrum of the baseband signal m(t)') ; 
xlabel('f (KHz)') ; ylabel('|M(f)|') ;  ylim([0 1.6]) ;
subplot(212)
plot(f/1000,abs(Sf),'r','LineWidth',2) ; grid ; title('Double-Sided Magnitude Spectrum of the baseband signal m(t)') ; hold on
plot(f_ds/1000,abs(Sf1),'b','LineWidth',1) ; grid ; title('Double-Sided Magnitude Spectrum of the sampled signal m(t)') ; ylim([0 1.8])
xlabel('f (KHz)') ; ylabel('|M1(f)|');

Note that I changed the colour and line-width of the plotting too to make the two curves more distinct.

HTH

[Mohamad]

I know from digital signal processing , when we sample the continuous time signal ( no aliasing ) , then the spectrum of the sampled signal will be a scaled and replicas of the spectrum of the original continous time signal , so I should have the original spectrum ( centererd at zero ) , it should appear without using the hold on in the script , i attach the | S1(f) | plot without hold on , i see only replicas but don't see the original spectrum , so any explainantion ?

[Bjorn Gustavsson]

What does "a scaled and replicas" mean?

If you call plot multiple times without hold on then you will only have the plot from the last call in your figure.

Now it is not clear to me what your problem is - this typically mean that you'll have to explain your problem more clearly.

[Mohamad]

I mean the following : I have xa(t) which is a continuous time signal it has a Fourier transfrm Xa(f) , then I sample xa(t) to get a sampled signal xs(t) then the Fourier transform for xs(t) is Xs(f) will consists of periodically repeated copies of Xa(f) , scaled by ( 1/T ) , These copies are shifted by integer multiples of the sampling frequency, and then superimposed to produce the periodic Fourier , given the equation below .

So I made sampling using downsample commnad , then fft , then why in my plot i didn't get the original spectrum ( centered at zero ) without using the hold on command ?

[Bjorn Gustavsson]

OK, now I see. That is how the discrete Fourier-transform works, in principle, with the different Nyquist zones. However, when you call matlab's (and as far as I know most other fft-implementations) with the Nfft-argument you get a zero-padding of your signal to that size. This changes the output from what you expect to what you get here.

Hopefully that clarifies the matter.

[Mohamad]

Yes , now I understand , I will not get the spectrum excatly as in the equation I attached before , thanks a lot .

Answer 2

[Paul]

Here is an illustration using a different signal; maybe it will provide some insight into your specific problem.

Define a continuous time signal

syms t real
f(t) = -(exp(-t*1i)*(exp(t*1i) - 1)^2)/t^2;

Plot f(t)

figure
fplot(abs(f(t)),[-30 30])
ylim([0 1])

We see that f(t) is quite small for abs(t) > 20. We'll use that fact later.

The Fourier transform of f(t) is

syms w real
F(w) = fourier(f(t),t,w)

F(w) = 

We see that F(w) is real. Plot it versus w.

figure
fplot(F(w),[-3 3])

We see that F(w) is a triangular pulse with cutoff frequency wc = 1;

Now, in order to not have aliasiing in the frequency domain after samplling, we have to choose a sampling period, Ts = 2*pi/w0, such that w0 > 2*wc. So let's choose w0 = 3*wc.

wc = 1;
w0 = 3*wc;
Ts = 2*pi/w0

Ts = 2.0944

We know that f(t) is of infinite duration, but based on the plot above we'll approximate f(t) by assuming f(t) = 0 for abs(t) > 10*Ts. With that in mind, the sampled signal is then approximated by the sum of the product of samples of f(t) and an impulse train. Note that f(0) = 1 which we have to define explicitly, otherwise we get a dvide by zero error.

f(t) = piecewise(t == 0, 1 ,f(t));
fs(t) = simplify(sum(f((-10:10)*Ts).*dirac(t-(-10:10)*Ts)));

Its Fourier transform is

Fs(w) = fourier(fs(t),t,w);

Convert it to a function for fast numerical computations

Fsfunc = matlabFunction(Fs(w));

Now we can evaluate and plot Fs as a function of frequency, and compare it to F(w)/Ts.

wvals = -10:.01:10;
figure
plot(wvals,abs(Fsfunc(wvals)))
hold on
fplot(F(w)/Ts)

As expected, the Fourier transform of the sampled signal is the sum of shifted replicas of the Fourier transform of the signal itself scaled by 1/Ts. I think the peaks of the blue curve a bit less than three and the little ripples around zero between the triangles result because we could only use a finite number of samples of f(t) in forming fs(t).

Now, can we get the blue curve using discrete-time processing? Yes, to good approximation. The continuous time Fourier transform (CTFT) of fs(t) is well-approximated by the discrete time Fourier transform (DTFT) of the samples of f(t). Because fs(t) is (approximately) finite duration, the sequence fs[n] formed from the samples of f(t) is also finite duration, and therefore we can use freqz() to compute its DTFT. However, freqz() assumes the first point in the sequence is at n = 0, whereas in our case it's at n = -10. So we have to use the time shifting property combined with freqz() to compute the DTFT. Note that one trip around the unit circle is equivalent to 2*pi/Ts = 3 rad/s, so we are actually computing the DTFT over several of its periods (recall that the DTFT is always 2*pi periodic).

figure
fsamples = double(f((-10:10)*Ts));
h = freqz(fsamples,1,wvals*Ts).*exp(1j*10*wvals*Ts);
plot(wvals,real(h)) % plote the real part because imag(h) is numerical error.

As expected, the DTFT of fs[n] is good approximation to the CTFT of fs(t).

If we want to use fft() to compute the DFT of the sampled sequence, it's very important to understand that, by definition, the DFT is non-zero only over one period around the unit circle. So we can use fft(), but we only get one period of the DTFT, which is the aprroximation of the CTFT of fs(t). Again, we have to deal with the time shift.

nfft = numel(fsamples);
wf1 = (0:nfft-1)/nfft*2*pi;
F1 = fft(fsamples,nfft).*exp(1j*10*wf1);
figure
plot(wf1/Ts,real(Fsfunc(wf1/Ts)))
hold on
stem(wf1/Ts,real(F1))

If we want to get the DFT samples more closely spaced, we zero pad

nfft = 64;
wf2 = (0:nfft-1)/nfft*2*pi;
F2 = fft(fsamples,nfft).*exp(1j*10*wf2);
figure
plot(wf2/Ts,real(Fsfunc(wf2/Ts)))
hold on
stem(wf2/Ts,real(F2))

Because the DFT is fundamentally limited to only cover one period of the DTFT, the DFT alone can't be used to illustrate the replication of the CTFT of f(t) that makes up the CTFT of fs(t) (it's a good approximation to samples of the CTFT of f(t)). Instead, we have to use repmat() on F2 (or F1), for example for four replicates.

figure
F3 = repmat(F2,1,4);
wf3 = (0:numel(F3)-1)/numel(F3)*4*2*pi/Ts;
stem(wf3,real(F3))

In summary, we can't use fft() to ilustrate how the CTFT of a sampled signal is the sum of shifted replicas of the CTFT of the signal itself. A bit of additional work is needed.

[Bjorn Gustavsson]

An even more direct illustration of this can be obtained by looking at this from the basics of the FT as a linear transform of a time-varying signal to a set of orthogonal basis-function that happens to be harmonically oscillating:

DFTMTX_full = dftmtx(numel(t));        % Full DFT-matrix, proper orthonormal
DFTMTX_down = DFTMTX_full(:,1:10:end); % The components corresponding to the Down-sampled signal, 
% no longer orthonormal due to much aliasing
fftDF = DFTMTX_down*sig_sampled(1:10:end)';
plot(f,(abs(fftshift(fftDF))))

This is not "that computationally efficient" - but a basic illustration of this as an effect of the aliasing of the no longer orthogonal harmonic basis-functions.

[Paul]

Hi @Bjorn Gustavsson

The orginal Question asked about the relationship between the CTFT of a signal and the CTFT after applying impulse modulation to that signal, all in the continuous-time domain. I only went to the disrete-time domain to show the relationship between the DTFT of a sequence of samples from the continuous-time signal and the CTFT of the coninuous-time, impulse-modulated signal, and that the DFT of the sampled sequence only covers a single period of that DTFT.

The code in the comment starts with a (finite-duration) sequence, sig_sampled, and compares the DFT of sig_sampled to the DFT of sig2, where sig2(n) = sig_sampled(n) for n = 1:10:end, and 0 otherwise. I don't understand how this analysis directly relates to continuous-time impulse modulation. I suppose it's similar in that the effect of impulse modulation in the continuous-time domain results in a signal that's zero between the impulses and this analysis is an analgous pulse modulation in the discrete-time domain.

[Bjorn Gustavsson]

@Paul: Your answer was a good an thorough explanation. My comment was intended as an additional illustration of the same effect, from another perspective, with the hope that it might illustrate how the aliasing-characteriestics of the DFT can be used for estimating the FT of band-limited signals from other frequency-bands than the base-band.

Also: downsample is a simple down-sampling similar to:

sig2 = sig_sampled(1:(n):end);

[Paul]

@Bjorn Gustavsson: Thanks for the clarification.

My understanding of all of the signals in this discussion is:

s(t) - continuous time signal

s1(t) - continuous time signal, impulse moudulated version of s(t) with modulation period T. The CTFT of s1(t) is constructed from a scaled sum of shifted replciates of the CTFT of s(t).

s1[n] - discrete time signal, s[n] = s(n*T)

s2[n] - downsampled sequence of s1[n], let M be the downsampling factor, s2 = s1(1:M:end) if s1 is finite duration. The DTFT of s2[n] can be constructed from a scaled sum of shifted and frequency-scaled replicates of the DTFT of s1[n]

s3[n] - upsampled sequence of s2[n], if s2[n] is finite duration, then s3(1:M:numel(s1)) = s2;

With these defintions and DFTMTX_down defined as above, this line

DFTMTX_down*s2(:)

yields frequency domain samples of the DTFT of s3[n].

[Bjorn Gustavsson]

@Paul: In the original code-snippet there are only 2 signals:

sig=2*cos(2*pi*400*t)+cos(2*pi*800*t)-3*sin(2*pi*1200*t); % original continuous time signal
sig_sampled = downsample(sig,fs/Fs);

That one obtains "aliasing-replicas" of the DFT for the downsampled sequence can be readily seen with this snippet:

DFTM = dftmtx(64);
subplot(2,2,1),
imagesc(real(DFTM)),
title('Full 64x64 DFT-matrix')
ylabel('Real part')
subplot(2,2,3),
imagesc(imag(DFTM))
ylabel('Imaginary part')
subplot(2,2,2),
imagesc(real(DFTM(:,1:4:end))),
title('Sownsampled in time 64x16 DFT-matrix')
subplot(2,2,4),
imagesc(imag(DFTM(:,1:4:end)))

[Mohamad]

Thanks a lot , great explainantion , the original signal is time limited , hence not band limited , hence i can not get excatly the spectra as the theoritcal eqaution ( i.e the spectra of the sampled signal ) .

[Paul]

@Bjorn Gustavsson: I fully agree with you that "aliasing-replicas" also occurs in discrete time and is an effect of downsampling. The only point that I was trying to make is that multiplication by DFTMTX_down captures the effect of upsampling the downsampled signal. For example:

s1 = 1:100;
DFTMTX_full = dftmtx(numel(s1));
M = 10;
s2 = downsample(s1,M);
s3 = upsample(s2,M);

DFTMTX_down = DFTMTX_full(:,1:M:end);

max(abs(fft(s3) - (DFTMTX_down*s2.').'))

ans = 2.9296e-14

In my view s3 results from pulse modulation of s1 in the discrete-time domain, which would be anlagous to impulse modulation in the continuous-time domain. The DTFT of s3 can be shown to be a the sum of shifted replicates of the DTFT of s1, which is analogous the sum of the shifted replicates of the CTFT in the continous-time domain.

% DTFT of S1
S1 = @(w) freqz(s1,1,w);

% DTFT of S3
S3 = @(w) freqz(s3,1,w);

w = 0:0.01:2*pi;

% create DTFT of s3 by sum of replicates of DTFT of s1
S3d = 0*S3(w);
for ii = 0:M-1
    S3d = S3d + S1(w - 2*pi*ii/M)/M;
end

plot(w,abs(S3(w)),'b',w,abs(S3d),'-o') % not sure why the solid line isn't blue?