%
% This script find the active filter of FM1-FM10
% and returns the total frequency response of the specific filter bank.
% 
% Modified by Masayuki Nakano, 2016-04-13
% Keiko Kokeyama, 2012-08-08
% originally by Tobin Fricke, 2010-04-19
%
% 
% - Checks which FMs are active
% - Calculate the total frequency response of the module
%
%
%
% input varargin :
% 1 the foton file name, string. prefer to use find_FilterFile function.
% (find_FilterFile('kamioka','K1','K1IMC',1144558290)...)
% 2 filter bank name, string.('K1:IMC-MCL_SERVO',...)
% 3 module name, string.('IMC_MCL_SERVO',...)
% 4 frequency vector, num vector
% 5 plot flag (true or false) to plot the total active module response
%
% output
% ModResp : frequency response of the bank module you specify (complex number)
%
%
% example
% frequency = logspace(-2, 4, 500);
% ModResp = ...
% getFilterResp(find_FilterFile('kamioka','K1','K1IMC',1144558290),...
%               'K1:IMC-MCL_SERVO','IMC_MCL_SERVO',logspace(0,3),true)


function ModuResp = getFilterResp(varargin)

% modified by M.Nakano
%ifo = varargin{1};
%subsystem = varargin{2};
%foton = varargin{3};
%bank_name = varargin{4};
%freq = varargin{5};
foton = varargin{1};
bank_name = varargin{2};
module_name = varargin{3};
freq = varargin{4};

sos = [];

%% Read the filter and state

filters = readFilterFile(foton);

%modified by M.Nakano
cmd1 = ['[a1, b1] = system(''z read ' bank_name '_SW1R '');'];
cmd2 = ['[a2, b2] = system(''z read ' bank_name '_SW2R '');'];
%cmd1 = ['[a1, b1] = system(''ezcaread -n ' [ifo] ':' [subsystem] '-' [bank_name] '_SW1R '');'];
%cmd2 = ['[a2, b2] = system(''ezcaread -n ' [ifo] ':' [subsystem] '-' [bank_name] '_SW2R '');'];

eval(cmd1)
eval(cmd2)

p = ezcaswitchreport(str2num(b1),str2num(b2));


for ii=1:length(p),
    newsos = filters.(module_name)(p(ii)).soscoef;
    sos = vertcat(sos, newsos);
end
fs = filters.(module_name)(1).fs;

% Choose the frequencies at which we want to know the response

% Call SOS2FREQRESP
ModuResp = sos2freqresp(sos, 2*pi*freq, fs);

%% Make a Bode plot

if varargin{5} == true
H = ModuResp;
subplot(2,1,1);
semilogx(freq, db(H));
xlim([min(freq) max(freq)]);
grid on;
ylabel('gain [dB]');
title(sprintf('Bode plot of all %s active filter modules', bank_name), 'interpreter','none');

subplot(2,1,2);
semilogx(freq, angle(H)*180/pi);
xlim([min(freq) max(freq)]);
set(gca, 'YTick', 45*(-4:4));
grid on;
ylabel('phase [degrees]');
xlabel('frequency [Hz]');
end