Spectra replica after sampling a continuous time signal
7 次查看(过去 30 天)
显示 更早的评论
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)|');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/919599/image.jpeg)
0 个评论
采纳的回答
Bjorn Gustavsson
2022-3-9
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
7 个评论
Bjorn Gustavsson
2022-3-10
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.
更多回答(1 个)
Paul
2022-3-12
编辑:Paul
2022-3-12
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)
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
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.
7 个评论
Paul
2022-3-17
编辑:Paul
2022-3-17
@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.').'))
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?
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Transforms 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!