# How to calculate step and impulse response at multiple sampling rates

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.

• I'm not sure why you chose this approach. Isn't the impulse response simply $h(t)=e^{-at}u(t)$? You can plot this at any sampling rate you like. And from $h(t)$ you can easily calculate the step response and plot it. Why do you want to use Matlab's conv() function? Jun 21, 2014 at 9:03
• I need to write some c++ code that approximates samples from the continuous time differential equation (described above). I need to make sure I'm implementing this correctly, so I was using matlab to prototype. I am also aware that FFT is usually preferred over direct convolution since it's $O(Nlog(N))$ vs $N^2$. Jun 23, 2014 at 17:14

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

• Decibels are not described by natural logarithm $\ln$, you must use the $\log_{10}$.
– jojek
Jun 23, 2014 at 17:27
• Impulse function was a posting error. Thank you for pointing it out. I wanted to write a procedure that accepted samples from x(nT) and output estimates of y(nT). I wanted T to be a parameter to the procedure. With your post and Phonon, I see it's not a trivial problem. Jun 23, 2014 at 17:30
• @jokek - My mistake. I was going to first talk about gain in DB's, but the 10 factor involved in the DB formula sometimes seems to confuse some people, so I decided to use $ln$ instead to make my point clearer, but I forgot to change the text. Thanks for pointing it out since it was a serious omission. Jun 24, 2014 at 14:11
• @Budric - This doesn't mean you should give up on the problem. Try to experiment with some simple functions and see what you get. Sometimes learning to interpret your results is more valuable that the results themselves, and you have some guidelines now on what to expect. See at least if what you get makes sense, and if not try to figure out why. Jun 24, 2014 at 14:24

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.

• Thank you, that makes more sense now. I think I need to investigate approximating differential equations in terms of difference equation to actually implement this. I've seen several mapping options from Laplace s domain to z-transform domain. Jun 23, 2014 at 17:19