H2ASQ - produce a stable linear filter that, applied to strain h,
produces AS_Q counts
[b,a] = h2asq(fname,nb,na,woff)
fname file containing frequency, in-phase and quadrature phase
components from, e.g., swept sine calibration, representing the
ratio of the excitation (e.g., ETMX) at that frequency to the AS_Q
excitation
nb, na filter order
woff invfreqz produces warnings of poorly conditioned matrices, though
it does appear to give acceptable results. Setting woff to a
non-zero value turns off the warning messages. This is here
because the user should know they are turning-off this warning!
b, a transfer function s.t. asq = filter(b,a,h);0001 function [b,a] = h2asq(fname,nb,na,woff) 0002 % H2ASQ - produce a stable linear filter that, applied to strain h, 0003 % produces AS_Q counts 0004 % 0005 % [b,a] = h2asq(fname,nb,na,woff) 0006 % 0007 % fname file containing frequency, in-phase and quadrature phase 0008 % components from, e.g., swept sine calibration, representing the 0009 % ratio of the excitation (e.g., ETMX) at that frequency to the AS_Q 0010 % excitation 0011 % nb, na filter order 0012 % woff invfreqz produces warnings of poorly conditioned matrices, though 0013 % it does appear to give acceptable results. Setting woff to a 0014 % non-zero value turns off the warning messages. This is here 0015 % because the user should know they are turning-off this warning! 0016 % 0017 % b, a transfer function s.t. asq = filter(b,a,h); 0018 0019 % $Id: h2asq.m,v 1.1 2003/09/24 18:59:13 lsf Exp $ 0020 % based on Landry's code to produce AS_Q to strain filters 0021 0022 % Deal with warning message 0023 if (nargin == 4) 0024 if (~isempty(woff)) 0025 if (woff ~= 0) 0026 warning off MATLAB:nearlySingularMatrix; 0027 end 0028 end 0029 end 0030 0031 % Read data file containing frequency, response information 0032 % Each line contains three numbers: 0033 % frequency (Hz) 0034 % in-phase component of ratio AS_Q/ETMX 0035 % quad-phase component of ratio AS_Q/ETMX 0036 [f, Re, Im]=textread(fname,'%f %f %f'); 0037 fs = 16384; % sampling frequency (Hz) 0038 w = 2*pi*f; % frequency in radians 0039 cr = Re + i*Im; % complex ratio AS_Q/ETMX; 0040 0041 0042 % find transfer function H(z) = B(z)/A(z) s.t. H(exp(2*pi*i*f/fs)) is equal 0043 % to Re + i*Im at each f 0044 % Emphasis to give each frequency in fitting the transfer function 0045 % Landry gives extra emphasis to frequencies in range [40,81] Hz. 0046 % LSF finds that we need to put more emphasis on frequencies in the 0047 % >400 Hz range to get a good fit in the gw band. 0048 wt = ones(size(cr)); 0049 %wt(find(40 < f & f < 81)) = 1.3; 0050 wt(find(400 < f)) = 1.5; 0051 % Note: this use of invfreqz guaranteed to produce stable filter 0052 iter = 10*length(f); % number of iterations for filter finding 0053 [b,a] = invfreqz(cr,w/fs,nb,na,wt,iter,[]); 0054 %[b,a] = invfreqs(cr,w,nb,na,wt,iter); 0055 0056 % compare 0057 [h,w] = freqz(b,a,w/fs); 0058 subplot(2,1,1); 0059 semilogx(f,20*log10(abs(cr)),f,20*log10(abs(h))); 0060 title(sprintf('[nb, na] = [%d, %d]',nb, na)); 0061 legend('Input data','Approx filter'); 0062 ylabel('dB'); 0063 grid; 0064 subplot(2,1,2); 0065 semilogx(f,abs(h)./abs(cr)); 0066 v = axis; 0067 v(3:4) = [1/2,2]; 0068 axis(v); 0069 xlabel('Frequency (Hz)'); 0070 ylabel('Ratio approx/actual'); 0071 grid; 0072 0073 % convert to zpg filter 0074 [zz,pp,kk]=tf2zp(b,a); 0075 0076 % add two poles to *excitation* so that excitation becomes strain 0077 % Because H is AS_Q/EXC, this is adding two zeros to H 0078 % Frequency of pendulum is so low compared to Nyquist that we don't need to 0079 % worry the difference between discrete, continuous filter 0080 zz = [zz(:); 2*pi*(0.076 + 0.76i)/fs; 2*pi*(0.076 - 0.76i)/fs]; 0081 0082 % now go back to transfer function form 0083 [b,a] = zp2tf(zz,pp,kk); 0084 0085 % we be done 0086 return