LIGO front-end digital decimation filter

README

This is the frequency response of the filter used before downsampling from 64 kHz to 2 kHz.

decim.m

%% Plot the digital decimation filter used in LIGO's controller.c
%
% Tobin Fricke
% 2012-03-26

% From controller.c:
%
% /* Coeffs for the 32x downsampling filter (2K system) per Brian Lantz May 5, 2009 */
% static double feCoeff32x[9] =
% 	{0.010581064947739,
%         -1.90444302586137,    0.91078434629894,   -1.96090276933603,    0.99931924465090,
%         -1.92390910024681,    0.93366146580083,   -1.84652529182276,    0.99866506867980};
% 

% The coefficients from the C code:

gain = 0.010581064947739;

sos = [ -1.90444302586137,    0.91078434629894,   -1.96090276933603,    0.99931924465090; ...
        -1.92390910024681,    0.93366146580083,   -1.84652529182276,    0.99866506867980];

% Put it into standard form
%
% The coefficients stored in the C code assume a form of the second order
% section where the leading coefficient in the numerator and the
% denominator are both 1.  Here we insert those two 1's into the matrix,
% and also swap the numerator and denominator.

rows = size(sos,1);
sos = [ ones(rows,1) sos(:,3:4) ones(rows,1) sos(:,1:2) ];

% Get the response
f = logspace(log10(10), log10(2^15),501);
H = gain * sos2freqresp(sos, 2*pi*f, 2^16);

a = subplot(2,1,1);
semilogx(f, db(H));
grid on;
ylabel('gain [dB]');

a(2) = subplot(2,1,2);
semilogx(f, unwrap(angle(H))*180/pi);
xlim([min(f) max(f)]);
set(gca, 'YTick', 45*(-8:4));
grid on;
ylabel('phase [degrees]');
xlabel('frequency [Hz]');
linkaxes(a,'x');

% subplot(2,1,2);
% semilogx(f, 1e6 * unwrap(angle(H))*180/pi  /360./f);
% xlim([min(f) max(f)]);
% grid on;
% ylabel('phase delay [us]');
% xlabel('frequency [Hz]');

title(a(1), 'feCoeff32x');
print -dpdf feCoeff32x.pdf
print -dpng feCoeff32x.png

feCoeff32x.pdf

feCoeff32x.png

G2FOO.txt

# FILTERS FOR ONLINE SYSTEM
#
# Computer generated file: DO NOT EDIT
#
# MODULES FeCoeff32x
#
# SAMPLING RATE 65536
#
################################################################################
### FeCoeff32x                                                               ###
################################################################################
# DESIGN   FeCoeff32x 0 sos(0.010581064947739,[-1.9609;0.999319;-1.90444;0.910784;-1.84653;0.998665;-1.92391;0.933661]) \
#                       
# DESIGN   FeCoeff32x 1 sos(0.010581064947739,[-1.9609;0.999319;-1.90444;0.910784]) \
#                       sos(1.000000000000000,[-1.84653;0.998665;-1.92391;0.933661]) \
#                       
###                                                                          ###
FeCoeff32x 0 21 2      0      0 AA              0.010581064947739  -1.9044399999999999   0.9107840000000000  -1.9609000000000001   0.9993190000000000
                                                                 -1.9239100000000000   0.9336610000000000  -1.8465300000000000   0.9986650000000000
FeCoeff32x 1 21 2      0      0 AA2lines        0.010581064947739  -1.9044399999999999   0.9107840000000000  -1.9609000000000001   0.9993190000000000
                                                                 -1.9239100000000000   0.9336610000000000  -1.8465300000000000   0.9986650000000000

sos2freqresp.m

function H = sos2freqresp(sos, omega, fs)
%SOS2FREQRESP  Frequency response of SOS filters.
%
%   H = SOS2FREQRESP(SOS,W,FS) computes the frequency response H of the 
%   matrix of second order section coefficients at the frequencies
%   specified by the vector W. These frequencies should be real and in
%   radians/second.   FS is the sample rate in samples per second.
%
%   Example usage:
%
%     filters = readFilterFile('L1FOO.txt');
%     sos = filters.('FOO_NOTCH')(1).soscoef;
%     fs = filters.('FOO_NOTCH')(1).fs;
%     f = logspace(log10(8), log10(12), 1001);
%
%     H = sos2freqresp(sos, 2*pi*f, fs); 
%
%     subplot(2,1,1);
%     semilogx(f, db(H));
%     subplot(2,1,2);
%     semilogx(f, angle(H)*180/pi);
%
% See also FREQRESP, ZP2SOS, SOS2ZP, SOS2SS, SS2SOS

% Author: Tobin Fricke, <tfricke@ligo.caltech.edu> 2010-04-19
% Louisiana State University and Agricultural and Mechanical College

% Validate the input
[rows, cols] = size(sos);
if cols ~= 6
    error('SOS matrix should have 6 columns');
end

% Do the computation
T = 1/fs;
s = 1i*omega;
z = exp(s * T);

% This should work no matter what the shape of omega.
H = ones(size(omega));

for ii=1:rows,
    b0 = sos(ii,1);
    b1 = sos(ii,2);
    b2 = sos(ii,3);
    a0 = sos(ii,4);
    a1 = sos(ii,5);
    a2 = sos(ii,6);
    
    H = H .* (b0 + b1./z + b2./z.^2) ./ (a0 + a1./z + a2./z.^2);
end

return;