What does the spectrum of a sampled signal look like if one uses a zero-order-hold?
There is some information in: Sampling theorem and Dirac comb, but I am wondering, how to apply this in practice.
Starting with an analytic spectrum I generated a time signal and recalculated the spectrum and it fits quite well. In a next step I used a moving average for the time signal and again calculated the spectrum. The moving average filter can be considered as a convolution with a rectangular window. For the Fouier Transform this is simply a multiplication and finally the analytical spectrum can be found by multiplying the original spectrum with the sinc function and its conjugate. The spectrum based on the data again fits quite well to the analytical spectrum.
But I am unable to find a similar solution for the signal sampled with the zero-order-hold. Some help would be highly appreciated :-)
Below the matlab code.
clear all; close all;clc
%% Configuration
T_max = 2^9; % [s]
f_max = 2^2; % [Hz]
fs = f_max*2^2;% [Hz] Nyquist–Shannon sampling theorem: fs>=f_max*2
SampleLength = 2^4; % [-]
%% Analytical Spectra
df = 1/T_max;
f = [df:df:f_max];
n_FFT = fs/df;
dt = 1/fs;
t = 0:dt:T_max-dt;
T = SampleLength*dt;
%-----------
S = 1./(1+f.^2);
%-----------
RECT = sinc(f*T);
X = 1./(1+1i*f);
S_MAV = X.*conj(X).*RECT.*conj(RECT);
%-----------
nk = 2^7; % T_max/T seems be a better number
for ik=-nk:1:nk
for jk=-nk:1:nk
X_i = 1./(1+1i*(f-ik/T));
X_j = 1./(1+1i*(f-jk/T));
S_ZOH_ij(jk+nk+1,:) = RECT.*conj(RECT).*X_i.*conj(X_j)*(1/T/4);
end
S_ZOH_i(ik+nk+1,:) = sum(S_ZOH_ij);
end
S_ZOH = sum(S_ZOH_i);
%% Time Signals
phi = rand(size(f))*2*pi; % random phase shift
A = sqrt(S/(1/2/df)); % Amplitudes
% time signal is a sum of all frequencies
x = (A*cos(2*pi*f'*t+repmat(phi',1,n_FFT)));
% Moving AVerage
x_MAV = filter(ones(1,SampleLength)/SampleLength,1,x);
% Zero Order Hold
x_ZOH = reshape(repmat(x(1:SampleLength:end),SampleLength,1),1,[]);
%% Spectra from data
nBlocks = 2^5;
nDataPerBlock = floor(n_FFT/nBlocks/2)*2; % should be even
vWindow = hann(nDataPerBlock);
%-----------
[S_D f_D] = pwelch(detrend(x,'constant'), vWindow,[],[],fs);
[S_D_MAV f_D] = pwelch(detrend(x_MAV,'constant'), vWindow,[],[],fs);
[S_D_ZOH f_D] = pwelch(detrend(x_ZOH,'constant'), vWindow,[],[],fs);
%% Spectra via FFT
x_FFT = fft(x,n_FFT);
X = [x_FFT(1) 2*x_FFT(2:n_FFT/2) x_FFT(n_FFT/2+1)]/n_FFT;
S_FFT = 1/2/df*X.*conj(X);
f_FFT = 0:df:1/2*fs;
%% Plots
figure('Name', 'Time')
hold on;grid on;box on;
plot(t,x,'.-');
plot(t,x_MAV ,'g');
plot(t,x_ZOH ,'o-r');
xlim([0 8])
xlabel('t [s]')
ylabel('x [-]')
legend('x','x_{MAV}','x_{ZOH}')
%-----------
figure('Name', 'Spectra')
hold on;grid on;box on;
plot(f_D,S_D,'LineWidth',3)
plot(f_D,S_D_MAV,'g','LineWidth',3)
plot(f_D,S_D_ZOH,'m','LineWidth',3)
plot(f,S,'c','LineWidth',1)
plot(f,S_MAV,'k','LineWidth',1)
plot(f,abs(S_ZOH),'r','LineWidth',1)
legend('S_{D}','S_{D,MAV}','S_{D,ZOH}','S','S_{MAV}','S_{ZOH}')
set(gca,'XScale','log')
set(gca,'YScale','log')
ylabel('S [1/Hz]')
xlabel('f [Hz]')
xlim([df*nBlocks f_max])
ylim([1e-4 1e1])
Update
I updated the code and we are getting closer (link). But there is still a problem with the scaling. Starting with the Fourier transform $$X(f)=\frac{1}{1+jf},$$ the auto spectrum is (omitting all scaling factors) $$S(f)=X(f)X^*(f)=\frac{1}{1+f^2},$$ where $X^*(f)$ is the conjugate, we get with Matt's comment using $\omega=2\pi f$ $$X_{ZOH}(f)=e^{-j\pi f T}sinc(fT)\frac{2\pi}{T}\sum_{k=-\infty}^{k=\infty} X(f-\frac{k}{T})$$ and finally $$S_{ZOH}(f)=X_{ZOH}(f)X_{ZOH}^*(f)=sinc^2(fT) \left(\frac{2\pi}{T}\right)^2 \left(\sum_{k=-\infty}^{k=\infty} X(f-\frac{k}{T})\right)\left(\sum_{k=-\infty}^{k=\infty} X^*(f-\frac{k}{T})\right).$$
For $k$ I used only a limited number, but the last two brackets are still computational intensive. It there any shortcut?
The remaining problem is the scaling: It seems to be in some way depending on the sampling frequency. Does anyone have an idea?
Update 2
By trial and error I found out, that the following Spectra is a good approximation with different sampling frequencies and hold times (link). $$S_{ZOH}(f)=X_{ZOH}(f)X_{ZOH}^*(f)=sinc^2(fT) \left(\frac{1}{4T}\right) \left(\sum_{k=-\infty}^{k=\infty} X(f-\frac{k}{T})\right)\left(\sum_{k=-\infty}^{k=\infty} X^*(f-\frac{k}{T})\right).$$ Does anyone have an explication for the $\left(\frac{1}{4T}\right)$ term?
Update 3
I added a part calculating the spectrum via FFT. At least for the unfiltered spectrum it coincide with the original spectrum. So, I assume, there is no scaling introduced by the pwelch algorithm.
The line A=sqrt(S/(1/2/df)) I used to get the amplitudes out of the spectrum to generate the time signal. I am not really sure, but it is related to the fact, that the spectrum is only one-sided and I use the power for one $df$: $$S(f)=\frac{A^2(f)}{2df}$$ Therefore, I don't think, it is the reason for the $\left(\frac{1}{4T}\right)$ term ...
