I am writing a software simulator to simulate a real physical acoustic interferometer. I now want to write the module to model the microphone amplifier(s). I'd like to be able to specify the frequency response of the amplifier both in terms of amplitude (gain) and phase. So I will feed in a sampled waveform and get a sampled waveform out that is modified according to the properties of the amplifier circuit.

I'm no expert in signal processing but happy to read the background maths. But to save me time going down dead end ideas could someone please tell me the standard method of doing this and what maths I have to read up on to be able to understand it?

  • $\begingroup$ What kind of amplifier is it? $\endgroup$ – Arnfinn Jan 2 '17 at 1:39
  • $\begingroup$ It's a general simulator. The system I'm modelling is a multi-channel sound recorder. I want to model the phase behaviour (as a function of frequency) through the whole system but the next step is the microphone amplifiers. I would like to input the frequency response in terms of amplitude and phase. I will then be able to study the effect of changing the amplifier properties with my simulator. $\endgroup$ – David Wallis Jan 2 '17 at 9:21
  • $\begingroup$ In practice it will be an op-amp or instrumentation amplifier. $\endgroup$ – David Wallis Jan 2 '17 at 9:27
  • $\begingroup$ Perhaps, after reading a little more, I can phrase the question a little better. I start with a waveform with 4096 samples. I have an amplifier with a known frequency response (amplitude and phase). How do I calculate the 4096 values that would come out of the amplifier? Perhaps this also boils down to: If I know the amplifier characteristics defined in terms of the frequency response (amplitude and phase), how do I define my filter (transfer function)? $\endgroup$ – David Wallis Jan 2 '17 at 11:51
  • $\begingroup$ Sounds like an amplifier with linear dynamics, so you can model it using an infinite impulse response (IIR) filter. Do you know how to implement an IIR filter realization? $\endgroup$ – Arnfinn Jan 2 '17 at 12:02

I apologize, the following is a bit rough, but I've put together an example in MATLAB illustrating all the different elements I think you need. That is:

  • a fist-order model for an amplifier (a common approximation for most op. amps.)
  • how to discretize it to find an IIR filter
  • one way to implement a realization of the IIR filter (straight forward state-space implementation)
  • how to generate a transfer-function estimate from measured data
  • one way (not the most robust, but perhaps the most obvious) to fit the transfer-function coefficients to the transfer-function estimate

To read up on the theory behind all this I would recommend the following text-books:

Chi-Tsong Chen, Linear System Theory and Design (3rd edition). Oxford University Press, 1999.

Petre Stoica and Randolph Moses, Spectral Analysis of Signals. Prentice Hall, 2004.

Lennart Ljung, System Identification: Theory for the User (2nd edition). Prentice Hall, 1999.

The last book is somewhat heavy, but concentrate on "output error" methods, and you should be very well covered for systems where you know the order of the model, or "subspace" methods if you don't. These methods are available in the MATLAB system identification toolbox.

%% Amplifier model
wb = 1e4; % Amplifier bandwidth
ga = 1; % Amplifier gain

% Modeling amplifier as having fist-order dynamics
% If using Octave run these in the command window:
%pkg install -forge control
%pkg load control
H = tf(ga,[1/wb 1]); % Continuous-time transfer function

%% Finding a realization (convert to state-space formulation)
% V/U = ga/(s/wb + 1)
% V*(s/wb + 1) = ga*U
% s*V = -wb*V + ga*wb*U
% dv/dt = -wb*v + ga*wb*u
% Standard form:
% dx/dt = A*x + B*u
% y = C*x + D*u
% with x = v, y = x
A = -wb; B = ga*wb; C = 1; D = 0;

%% Finding a discretization:
Ts = 1/1e5; % sampling time
% Using Euler integration
Ad1 = 1 + Ts*A; Bd1 = Ts*B; Cd1 = C; Dd1 = D;
% Using exact discretization (use when system has analog-digital converters)
Ad2 = exp(A*Ts); Bd2 = (Ad2 - 1)/A*B; Cd2 = C; Dd2 = D;

%% Computing a step-response to demonstrate one way to implement an IIR filter
u = [zeros(1,100) ones(1,200)]; % input
x = zeros(2,numel(u)+1);
y = zeros(2,numel(u));

t = 0:Ts:(numel(u)-1)*Ts;

% Comupting filter outputs
for k = 1:numel(u)
    x(1,k+1) = Ad1*x(1,k) + Bd1*u(k); % response of Euler integration scheme
    y(1,k) = Cd1*x(1,k) + Dd1*u(k);

    x(2,k+1) = Ad1*x(2,k) + Bd1*u(k); % response of exact discretization
    y(2,k) = Cd1*x(2,k) + Dd1*u(k);


%% Computing the response to noise for identification of parameters
% The frequency response of a system is in general most reliably found when
% exciting the system with white noise. A pseudo-random binary sequence is
% theoretically optimal when you don't know anything about your system, but
% can in practice cause problems due to the switched nature. A Gaussian
% noise source is usually a good substitute.
un = randn(1,1e4); % input
xn = zeros(1,numel(un)+1);
yn = zeros(1,numel(un));

% Comupting filter outputs
for k = 1:numel(un)
    xn(k+1) = Ad2*xn(k) + Bd2*un(k);
    yn(k) = Cd2*xn(k) + Dd2*un(k);
% Remove transient as you want to look at stationary response
un_ = un(201:end);
yn_ = yn(201:end);

% Compute power spectral densitites and cross-power spectral densitites
% using a simplified version of the Welch method
N = length(yn_); % length of original sequence
L = 5; % number of PSD averages
M = floor(N/L); % length of sequence segments
W = 0:(2*pi)/M:2*pi-(2*pi)/M; % normalized PSD frequencies

WIN = hanning(M); WIN = WIN(:).'; % optional windowing (if the segments are short)
Pwin = 1/M*sum(abs(WIN).^2); % window "power" correction (for Welch method)

% PSD for input sequence
Puu = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';
    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puu = Puu + (Uft.*conj(Uft));
Puu = 1/(L*Pwin)*Puu;

% Cross-PSD between input and output sequence
Puy = zeros(1,M);
for k = 1:L
    un_segment = un_((k-1)*M+1:k*M);
    un_segment = un_segment(:).';

    yn_segment = yn_((k-1)*M+1:k*M);
    yn_segment = yn_segment(:).';

    yn_windowed = yn_segment.*WIN;
    Yft = 1/sqrt(2*pi*M)*fft(yn_windowed);

    un_windowed = un_segment.*WIN;
    Uft = 1/sqrt(2*pi*M)*fft(un_windowed);

    Puy = Puy + (Yft.*conj(Uft));
Puy = 1/(L*Pwin)*Puy;

% Transfer-function estimate
Huy = Puy./Puu;


%% Fitting a model to the output data

% finding the z-transform of discrete-time model (Euler integration)
% using the standard expression
z = tf('z',Ts);
Hd1 = Cd1*(z - Ad1)\Bd1 + Dd1; 

% finding the z-transform of discrete-time model (exact discretization)
% using an alternative method
b0 = Bd2*Cd2;
a0 = -Ad2;
a1 = 1;
Hd2 = tf(b0,[a1 a0],Ts);

% compare the continuous and discrete-time models

zv = exp(1i*W); % starred transform (between z- and s-domains)

% H = b0/(a1*z + a0)
% H(W) = b0/(a1*exp(i*W) + a0)
% a1 = 1; % in order to have unique coefficients for a0 and b0
% H(W) = b0/(exp(i*W) + a0)
% (exp(i*W) + a0)*H(W) = b0
% exp(i*W)*H(W) = -a0*H(W) + b0
% exp(i*Ts*w)*H(W) = [-H(W); 1]*[a0, b0]
% zeta = phi*theta;

zeta = zv(:).*Huy(:);
phi = [-Huy(:), ones(size(Huy(:)))];

theta_true = [a0; b0]
theta_est = phi\zeta

The amplifier has a frequency response in the regular sense of the term only when it is operating in the linear region, for inputs of known and constant bandwidth. You may be able to piece it together from the spec sheet. The gain will depend on the gain-bandwidth product and the amplifier configuration; the phase may be stated in terms of time delay.

If you want to model the amplifier non-linearities, and depending on the accuracy you want in your model, this can be easy or very hard. I'll describe the basics of a simple approach, which may help to get you started.

An amplifier is a non-linear system. Many memoryless non-linear systems can be modeled as a power series. If $x(t)$ is the input and $y(t)$ is the output, then $$y(t)=\sum_{k=0}^N a_k x^k(t).$$ This model allows you to study many non-linear effects such as intermodulation, saturation, and AM-AM & AM-PM distortion.

The problem is deciding on the values of $N$ and the coefficients $a_k$. Most amplifier models I've seen use $N=3$, since this allows modeling of 3rd-order intermodulation products. The values of the coefficients may be harder to obtain; you may have to measure them experimentally.

  • $\begingroup$ So the output is a function of a polynomial. I hadn't considered a continuous model; I was considering that my solution would be a digital filter, with the phase modelled with complex quantities. I would presumably need to determine the coefficients by fitting to some data. What would you suggest I use as an input to the amplifier to get the best estimate of the coefficients? An impulse? A square wave would give me a broad spectrum signal to work. I was also imagining that the input that I would use to characterise the amplifier model would be the frequency response (gain and phase). $\endgroup$ – David Wallis Dec 30 '16 at 10:45
  • $\begingroup$ @DavidWallis An ampifier can't be modeled as a (LTI) digital filter. You can still discretize the power series above, but you need to make sure the sampling frequency is large enough to capture all the frequencies in the output without aliasing. I'm afraid that how to find the power series coefficients for a particular amplifier is beyond my grasp of the subject. $\endgroup$ – MBaz Dec 30 '16 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.