Spectrum of a Signal Going Through a Zero Order Hold Sampling System

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 ...

• I've changed the frequency $\omega$ in my answer to $f$, which has some consequences on the scaling. Hope it helps. By the way, where does your spectrum $X(f)=\frac{1}{1+jf}$ come from? Which time signal in your program does it correspond to? (sorry, probably I just haven't taken enough time to understand your code). May 20 '13 at 11:29

The sampled time domain signal with a zero-order hold is a weighted pulse train convolved with the impulse response of the zero-order hold:

$$\sum_n x(nT)\delta(t-nT) * \text{rect}(t-T/2,T/2)$$

We have the following Fourier transform pairs

$$\sum_n x(nT)\delta(t-nT) \Longleftrightarrow \frac{1}{T}\sum_k X(f-\frac{k}{T})\\ \text{rect}(t-T/2,T/2)\Longleftrightarrow e^{-j\pi fT}T\text{sinc}(fT)$$ with $\text{sinc}(x)=\sin(\pi x)/(\pi x)$.

Multiplying the two transforms, we get for the Fourier transform of the sampled signal with zero-order hold:

$$e^{-j\pi fT}\text{sinc}(fT)\sum_k X(f-\frac{k}{T})$$

The frequency response of the zero-order hold attenuates the shifted versions of the spectrum $X(f)$ somewhat but does not really suppress them well.

• Thanks Matt! Please find the update and the remaining problem above. May 19 '13 at 14:07
• I've changed the frequency $\omega$ in my answer to $f$, which has some consequences on the scaling. Hope it helps. May 20 '13 at 11:16
• Thanks Matt! Unfortunately I am still wondering, where the $\left(\frac{1}{4T}\right)$ term comes from. Based on your equations it should be only the $sinc^2$ and the two sums. Or did I made some mistakes? The spectrum $S(f)=\frac{1}{1+f^2}$ is just an example. The spectrum I am working with it more complicated. The time signal in the example is generated out of the spectrum. May 20 '13 at 16:18
• You should also check whether the function $\tt{pwelch}$ introduces some scaling. Since it is a power spectrum estimation I wouldn't be surprised if there is a $1/T$ factor which differs from your energy computation of deterministic signals. May 20 '13 at 18:13
• Doesn't the factor $1/4$ come from the line $\tt{A = sqrt(S/(1/2/df));}$? May 20 '13 at 18:18