Frequency and impulse response

Frequency response functions

Consider a dynamic system described by the input-output differential equation—with variable \(y\) representing the output, dependent variable time \(t\), variable \(u\) representing the input, constant coefficients \(a_i,b_j\), order \(n\), and \(m \le n\) for \(n\in\mathbb{N}_0\)—as: \[\begin{aligned} {8} \label{eq:ioode_frf} \frac{d^n y}{d t^n} & {}+{} & a_{n-1} \frac{d^{n-1} y}{d t^{n-1}} & {}+{} & \cdots & {}+{} & a_1 \frac{d y}{d t} & {}+{} && a_0 y =\nonumber \\ b_m \frac{d^m u}{d t^m} & {}+{} & b_{m-1} \frac{d^{m-1} u}{d t^{m-1}} & {}+{} & \cdots & {}+{} & b_1 \frac{d u}{d t} & {}+{} && b_0 u. \end{aligned}\]

The Fourier transform \(\mathcal{F}\) of yields something interesting (assuming zero initial conditions): \[\begin{aligned} {8} \mathcal{F}\left( \frac{d^n y}{d t^n}\right. & {}+{} & a_{n-1} \frac{d^{n-1} y}{d t^{n-1}} & {}+{} & \cdots & {}+{} & a_1 \frac{d y}{d t} & {}+{} && \left. a_0 y \vphantom{\frac{dx^2}{dt^2}} \right) =\nonumber \\ \mathcal{F}\left( b_m \frac{d^m u}{d t^m}\right. & {}+{} & b_{m-1} \frac{d^{m-1} u}{d t^{m-1}} & {}+{} & \cdots & {}+{} & b_1 \frac{d u}{d t} & {}+{} && \left. b_0 u \vphantom{\frac{dx^2}{dt^2}} \right) \quad \Rightarrow \\ \mathcal{F}\left( \frac{d^n y}{d t^n} \right) & {}+{} & a_{n-1} \mathcal{F}\left( \frac{d^{n-1} y}{d t^{n-1}} \right) & {}+{} & \cdots & {}+{} & a_1 \mathcal{F}\left( \frac{d y}{d t} \right) & {}+{} && a_0 \mathcal{F}\left( y \right) =\nonumber \\ b_m \mathcal{F}\left( \frac{d^m u}{d t^m} \right) & {}+{} & b_{m-1} \mathcal{F}\left( \frac{d^{m-1} u}{d t^{m-1}} \right) & {}+{} & \cdots & {}+{} & b_1 \mathcal{F}\left( \frac{d u}{d t} \right) & {}+{} && b_0 \mathcal{F}\left( u \right) \quad \Rightarrow \\ (j\omega)^n Y & {}+{} & a_{n-1} (j\omega)^{n-1} Y & {}+{} & \cdots & {}+{} & a_1 (j\omega) Y & {}+{} && a_0 Y =\nonumber \\ b_m (j\omega)^m U & {}+{} & b_{m-1} (j\omega)^{m-1} U & {}+{} & \cdots & {}+{} & b_1 (j\omega) U & {}+{} && b_0 U. \end{aligned}\] Solving for \(Y\), \[\begin{aligned} Y &= \frac{ b_m (j\omega)^m + b_{m-1} (j\omega)^{m-1} + \cdots + b_1 (j\omega) + b_0 }{ (j\omega)^n + a_{n-1} (j\omega)^{n-1} + \cdots + a_1 (j\omega) + a_0 } U.\end{aligned}\] The inverse Fourier transform \(\mathcal{F}^{-1}\) of \(Y\) is the forced response. However, this is not our primary concern; rather, we are interested to solve for the frequency response function \(H(j\omega)\) as the ratio of the output transform \(Y\) to the input transform \(U\), i.e.¹ $$\begin{align} \label{eq:frf_definition} H(j\omega) &\equiv \frac{Y(\omega)}{U(\omega)} \\ &= \frac{ b_m (j\omega)^m + b_{m-1} (j\omega)^{m-1} + \cdots + b_1 (j\omega) + b_0 }{ (j\omega)^n + a_{n-1} (j\omega)^{n-1} + \cdots + a_1 (j\omega) + a_0 }. \end{align}$$

Note that a frequency response function can be converted to a transfer function via the substitution \(j\omega \mapsto s\) and, conversely, a transfer function can be converted to a frequency response function² via the substitution \(s \mapsto j\omega\), as in \[\begin{aligned} H(j\omega) = H(s)|_{s\rightarrow j\omega}.\end{aligned}\] It is often easiest to first derive a transfer function—using any of the methods described, previously—then convert this to a frequency response function.

Frequency response

From above, we can solve for the output response \(y\) from the frequency response function by taking the inverse Fourier transform: \[\begin{aligned} y(t) &= \mathcal{F}^{-1}Y(\omega). \end{aligned}\] From the definition of the frequency response function , \[\begin{aligned} y(t) &= \mathcal{F}^{-1}(H(j\omega) U(\omega)). \end{aligned}\]

The convolution theorem states that, for two functions of time \(h\) and \(u\), $$\begin{align} \mathcal{F}(h*u) &= \mathcal{F}(h) \mathcal{F}(u)\\ &= H(j\omega) U(\omega), \label{eq:conv-01} \end{align}$$ where the convolution operator $*$ is defined by $$\begin{align} \label{eq:conv-02} (h*u)(t) &\equiv \int_{-\infty}^\infty h(\tau) u(t-\tau)\, d\tau. \end{align}$$ Therefore, $$\begin{align} y(t) &= \mathcal{F}^{-1}(H(j\omega) U(\omega)) \\ &= (h*u)(t) \tag{from \eqref{eq:conv-01}} \\ &= \int_{-\infty}^\infty h(\tau) u(t-\tau)\, d\tau. \tag{from \eqref{eq:conv-02}} \end{align}$$ This is the frequency response in terms of all time-domain functions.

Impulse response

The frequency response result includes an interesting object: \(h(t)\). What is the physical significance of \(h\), other than its definition, as the inverse Fourier transform of \(H(j\omega)\)?

Consider the singularity input \(u(t) = \delta(t)\), an impulse. The frequency response is \[\begin{aligned} y(t) &= \int_{-\infty}^\infty h(\tau) \delta(t-\tau)\, d\tau. \end{aligned}\] The so-called sifting property of \(\delta\) yields \[\begin{aligned} y(t) &= h(t). \end{aligned}\] That is, \(h\) is the impulse response.

A very interesting aspect of this result is that \[\begin{aligned} H(j\omega) = \mathcal{F}(h). \end{aligned}\] That is, the Fourier transform of the impulse response is the frequency response function. A way to estimate, via measurement, the frequency response function (and transfer function) of a system is to input an impulse, measure and fit the response, then Fourier transform it. Of course, putting in an actual impulse and fitting the response, perfectly are impossible; however, estimates using approximations remain useful.

It is worth noting that frequency response/transfer function estimation is a significant topic of study, and many techniques exist. Another method is described in .

sys = tf([1,20],[1,4,20])

sys =
 
      s + 20
  --------------
  s^2 + 4 s + 20
 
Continuous-time transfer function.

poles = pole(sys)

poles =

  -2.0000 + 4.0000i
  -2.0000 - 4.0000i

abs(poles(1))

ans =

    4.4721

fs = 1000; % Hz .. sampling frequency
N = 2^12;
t_a = 0:1/fs:(N-1)/fs;
h_a = impulse(sys,t_a);

noise = 0.01*randn(N,1);
h_noisy = h_a + noise;

figure
plot(...
  t_a,h_noisy, ...
  'linewidth',1.5 ...
)
xlabel('time (s)')
ylabel('impulse response')

H = fft(h_noisy);

H_mag = abs(H/fs); % note the scaling
H_mag = H_mag(1:N/2+1); % first half, only

H_pha = angle(H); % note the scaling
H_pha = H_pha(1:N/2+1); % first half, only

f = fs*(0:(N/2))/N;

figure
semilogx(...
  2*pi*f,H_mag, ...
  'linewidth',1.5 ...
)
xlabel('frequency (rad/s)')
ylabel('|H(j\omega)|')

figure
semilogx(...
  2*pi*f,180/pi*H_pha, ...
  'linewidth',1.5 ...
)
xlabel('frequency (rad/s)')
ylabel('\angle H(j\omega) (deg)')