How to calculate step and impulse response at multiple sampling rates

Question

I'm struggling with a basic concept. Suppose I have a system $y(t)=e^{-at}u(t)\ast x(t)$ and I want to plot impulse response and step response using Matlab conv() function at various sampling rates. I'm having trouble understanding what gain I should apply and why. I tried:

a = 0.05;
cutoffTime = 250;
T = 1/samplingRates(i);  %assume this is defined somewhere [1, 10, 20]
timeRange = 0:T:cutoffTime;
step = ones(1, length(timeRange));
impulse = zeros(1, length(timeRange));
impulse(1) = 1;
exponential = exp(-a*timeRange);
stepResponse = conv(exponential, step, 'full');  
impulseResponse = conv(exponential, step, 'full');

I understand that mathematically sampling is modeled as multiplication by a pulse train, which introduces gain of $1/T$ to the discrete system. So the reconstruction filter needs gain of T. In the above code, for different sampling rates (say 1, 10, 20 samples/second) if I don't apply the gain of T to step response, I get different sized step response. However if I apply the gain T to impulse response, I don't get the impulse response I expect for all 3 sampling rates. Can anyone clarify where I'm going wrong in my implementation?

Thank you.

Answer 1

To talk about gain you have to be able to compare the RMS values of your discrete and continuous functions, and the discrete function will only have an RMS value if you choose a way to connect the points of your discrete signal in $\mathbb{R}^2$. We can use rectangles since most DACs maintain the voltage constant in between samples, and in the limit when the sample rate tends to infinity it converges to the function.

Then the RMS value of your set of points $\{{x_1,x_2, ...,x_n}\}$ would be the square root of the summation of your $x_i^2$ times the length of the base of the rectangle $1/f_s$ divided by the total length of the interval $n/f_s$:

$x_{_{RMS}}=\sqrt{\frac{1}{n}\sum_{i=0}^{n}x_i^2}$

On the other hand, the RMS value for your continuous signal, calculating the area of the squared function divided by the length of the interval, and considering the beginning of the interval is $t = 0$ is given by

$f_{_{RMS}}=\sqrt{\frac{1}{t}\int_{0}^{t}[f(t)]^2dt}$

So you would have to find a normalization constant that makes the gain 0. It means:

$ 0=\ln{(\frac{\alpha x_{_{RMS}}}{f_{_{RMS}}})} \Rightarrow \frac{\alpha x_{_{RMS}}} {f_{_{RMS}}}=1 \Rightarrow \frac{f_{_{RMS}}}{ x_{_{RMS}}} =\alpha$

where $\alpha$ is your normalization constant, and as you can see it would depend on the function involved. However, it is important to notice that for the same time interval, $x_{_{RMS}}$ may change if the sampling rate changes, so the normalization value (if it exists, since f might not be integrable in the chosen interval) may depend on the sample rate, tending to be closer to 1 as the sample rate increases. The normalization constant for most functions will in fact depend on the chosen interval, so thinking in terms of a normalization constant is just an approximation. In particular calculating the gain of a filter is not a simple problem, and in many cases the gain depends on the frequency, so you choose a particular frequency to estimate the gain of the filter.

You say that the responses are not what you expect, but maybe providing more details on what you expect and what are you obtaining would allow us to provide you with a better answer.

Also in your code, you are using the step function for both the impulse response and step response, and your impulse response has only zeroes. You have to set one value to 1.

Also keep in mind that when analyzing impulse and step responses of a filter the way you are doing it, it is a common practice to use sample period as the time unit and not seconds, and the units for the frequency response would then be in terms of sampling frequency so you have a more general idea of the response of the filter.

Hope this helps. DP

Answer 2

This is not the full answer you want, but the problem is that you can't permute sampling and convolution. Let's say that the $D(\cdot)$ operator takes a continuous signal and returns its sampled version. What you're doing is

$$y'(n) = D\left(e^{-at}u(t)\right) * D\left(x(t)\right)$$

What you want to get is

$$D\left(y(t)\right) = D\left(e^{-at}u(t) * x(t)\right) \neq y'(n).$$

It is in general impossible to do what you're trying to do, but I'm sure that good approximations exist. I'm not an expert on what they are or how to use them.