# Fast algorithm for the MATLAB's $\tt impz$ function? [closed]

• Does anybody here know what is the code algorithm behind the impz function?
• I know it uses the recursive algorithm but how is such implemented in impz that it can get the impulse response of very high $N$ and $M$ orders (numerator and denominator coefficients)?

I can make a script that gets the impulse response of a second order LCCDE (linear constant coefficient differential equation) but I want to know a GENERAL code that can get the impulse response of any LCCDEs regardless of their order.

• Is this a Matlab function? Perhaps give some background for those who dont' know the particularities of the function?
– MBaz
Dec 2, 2016 at 0:02
• yes. The impz function generates the impulse response of ANY n order of a filter given its numerator and denominator. I want to know how is this to be implemented in a self made code. All I can do is just the impulse response of a 2nd order filter since it's just manual indexing. Dec 2, 2016 at 0:08
• I'm voting to close this question as off-topic because, again, this is a question easily answerable by the "Algorithm" section in the official Matlab docs. Dec 2, 2016 at 11:36
• @MarcusMüller the question is not something about matlab documentation but about a fast implementation of an esixting function. It may not be the most exciting question but I belive it makes sense to ask it. As you already know matlab is a valid and legal tag here... Dec 2, 2016 at 11:49
• @Fat32 "Does anybody here know what is the code algorithm behind the impz function?" certainly is a question answered in the docs. The question how to do it fast is only in the title, and I interpret the question in a manner that OP asks how the impz function is so fast Dec 2, 2016 at 13:23

the impz function of matlab, simply computes N samples of the impulse response of a given LTI system from its LCCDE description coefficients a nd b. As far as I know its efficiency is a secondary issue, since computing an impulse response in real time is not a frequently encountered task.

One equivalent code is simply this:

  % given a,b vectors and number of samples as N
x = [1 zeros(1,N-1)];
h = filter(b,a,x);


I guess you know how the filter function is implemented and what its output means. In the simplest terms, the filter() matlab function computes the output vector y for an input vector x of an LTI system: $\sum_{k=0}^{k=N}a_k y[n-k] = \sum_{k=0}^{k=M} b_k x[n-k]$ , by the recursive approach:

$$y[n] = - \sum_{k=0}^{k=N} \frac{a_k}{a_0} y[n-k] + \sum_{k=0}^{k=M} \frac{b_k}{a_0} x[n-k]$$

Where the initial conditions of the system is assumed to be zero, inline with impulse response definition and LTI system with initial rest conditions.

Note that in the impulse response computation, the signal $x[n]$ has only a single nonzero entry at $n=0$ and the rest of the computation is a recursion on the signal $y[n]$ alone.

• Thank you! Can I use freqz to find the impulse response of a given filter? I mean from freqz what can I do to retreive the impulse response? Dec 2, 2016 at 10:21
• This is "signal processing" stack exchange, not "teach me how to read matlab documentation and examples", sorry. Dec 2, 2016 at 11:39
• Yes you can use freqz() to find out the impulse response from given b,a vectors. It's pretty straightforward, as freqz() produces the first half of the Frequency Response $H(e^{j\omega})$ from $\omega =0$ to $\omega = \pi$ , you should first convert that into the full spectrum from $\omega =0$ to $\omega = 2\pi$ and perform an inverse FFT of appropriate length to get the proper impulse response $h[n]$ . As Marcus states such details can get out of the theory of signal processing however. Nevertheless matlab is a legal tag here and some questions are sure to be allowed Dec 2, 2016 at 11:46
• Thank you!! How many n points should I use in y= freqz(b,a,n,'whole')? And what if I want to get 100 samples of the impulse response. Should I use fft(y, 100)? Sorry if I question too much. i tried doing y = freqz(b,a,100,'whole') to divide the WHOLE circle by 100. Then do ifft(y,100) to get 100 h(n) samples but the results is different when I used impz(b.a). Dec 2, 2016 at 12:38
• you shall consider issues involved in FFT and inverse FFT usage and the signal lengths... Dec 2, 2016 at 15:18