Let's run our exercise on the empirical estimation of Transfer Functions on the Matlab terminal.

The idea of the exercise is the following:

  1. simulate some white noise x(t)
  2. build a low-pass filter F from a pole-zero model
  3. pass the input noise x(t) through the filter and add some more noise yn(t) at the output so to have y = F*x(t) + yn(t)
  4. evaluate and plot the transfer function x -> y
In a flow diagram, the representation is as follows:

Dataflow for the 1St example of ao/tfe

The command-line sequence is the following:

%% General definitions
nsecs = 10000;
fs = 1;
%% Input noise
x = ao(plist('waveform', 'noise', 'sigma', 1, 'fs', fs,  'nsecs', nsecs, 'yunits', 'V'))

%% Filter 
pzm = pzmodel(plist('name', 'lp_filter', 'poles', [pz(0.05)], ...
              'gain', 10, 'iunits', 'V', 'ounits', 'A'))
lp_filter = miir(plist('pzmodel', pzm, 'fs', fs))
xf = simplifyYunits(filter(x, lp_filter))

%% Output noise
yn = ao(plist('waveform','noise', 'sigma', 3, 'fs', fs, 'nsecs', xf.nsecs, 'yunits', 'A'))
y = xf + yn

%% Plotting input and output noise
xx = psd(x, plist('scale',...
                  'ASD',...
                  'nfft', 1000))
yy = psd(y, plist('scale','ASD', ...
                  'nfft', 1000))
iplot(xx, yy, plist('Arrangement', 'subplots', 'YRanges', {[1e-1 1e1], [1e-2 1e2]}));

Ohmic admittance between y and x

Now we can proceed with the call to the ao/tfe method. The parameter list is very similaer to the one employed for the other spectral estimators:

Key Value Description

NFFT

1000

The number of samples defining the lenght of the window to apply

WIN

'BH92'

Or a different one, if you want.

OLAP

-1

Overlap will be chosen based on the window properties

ORDER

0

Segment-wise detrending up to order 0

The command line is the following:

%% Estimate the x->y transfer function
tf = tfe(x, y, plist('nfft', 1000, 'win', 'BH92'), 'olap', -1, 'order', 0);
tfxy = index(tf, 1, 2);

Please note the usage of the
ao/index
method to access the elements of an aos matrix [2x2] in this case, without breaking the history.

We also would like to evaluate the expected transfer function, which is obviously the filter tranfer function, or response. This is calculated by means of the
miir/resp
method. A detailed description of digital filtering is available in the User Manual dedicated section and will be touched upon in this topic; here let's just use the simplest form, where the needed parameter is a list of the frequency to evaluate the response at:

Key Value Description

F

tfxy.x

a vector of frequency values or an AO whereby the x-axis is taken for the frequency values

So we can just pass the x field of the fsdata ao containing the transfer function estimate. The command line is the following:

%% Evaluate the expected x->y transfer function
rf = resp(lp_filter, plist('f', tfxy.x))
  

Eventually let's look at the results:

%% Plotting estimated and expected transfer functions
iplot(tfxy, rf, plist('colors',{[1 0 0],[0 0 0]},'YRanges', {[1e-1 1e2], [-200 200]}))

TFE of x into y