% Calc_FRF.m
%
% Calculates and displays the frequency response function from the 
% recorded signals: input (left) and output (right).

% Recorded measurement wavefiles
fNameP = 'gwn3vp';
fNameI = 'gwn3vi';
tText = 'noise';

% Parameter settings and initialisation
R_series = 0.5;
nW = 1024;
nA = nW / 2;
nZ = 8192;
DataP = zeros(nZ,2);
DataI = zeros(nZ,2);
win = hanning(nW);
fMax = 6;% maximum frequency displayed on the x-axis (kHz)

% Read in the stereo signals
[yP,fs] = wavread([fNameP '.wav']);
nTP = length(yP);
yP = [zeros(nW/2,2); yP; zeros(nW/2,2)];
[yI,fs] = wavread([fNameI '.wav']);
nTI = length(yI);
yI = [zeros(nW/2,2); yI; zeros(nW/2,2)];

% Loop through sections of the pressure and voltage waveforms
nX = ceil(nTP/nA);
for iX = 1:nX,
  % Set up a vector that includes some zero padding
  data = zeros(nZ,2);
  % Copy the current section of the waveform
  data(1:nW,:) = yP([1:nW] + (iX-1) * nA,:) .* [win win];
  % Add its contribution to the power spectrum
  DataP(:,1) = DataP(:,1) + abs(fft(data(:,1))).^2;
  DataP(:,2) = DataP(:,2) + abs(fft(data(:,2))).^2;
end% for
% Divide by the number of sections and the number of bins
DataP = DataP ./ (nX * nW);

% Loop through sections of current and voltage waveforms
nX = ceil(nTI/nA);
for iX = 1:nX,
  % Set up a vector that includes some zero padding
  data = zeros(nZ,2);
  % Copy the current section of the waveform
  data(1:nW,:) = yI([1:nW] + (iX-1) * nA,:) .* [win win];
  % Add its contribution to the power spectrum
  DataI(:,1) = DataI(:,1) + abs(fft(data(:,1))).^2;
  DataI(:,2) = DataI(:,2) + abs(fft(data(:,2))).^2;
end% for
% Divide by the number of sections and the number of bins
DataI = DataI ./ (nX * nW);

% Calculate the ratio of the output to the input
GainP = DataP(:,2) ./ DataP(:,1);
GainI = DataI(:,2) ./ (R_series.^2 * DataI(:,1));
GainA = GainP ./ GainI;

% Set up the scale for the frequency axis
freq = [0:nZ-1]' * fs / nZ;

% Plot the resulting gain in decibels versus frequency
figure(1);
subplot(211);
plot(freq/1000, 10*log10(abs([GainP])),'b-');
hold on;
plot(freq/1000, 10*log10(abs([GainI])),'g--');
hold off;
grid on;
xlabel('Frequency (kHz)');
ylabel('Gain (dB)');
title(['Frequency response functions per volt, with ' tText ' signal']);
axis([0 fMax -30 20]);
legend('P','I');

subplot(212);
plot(freq/1000, 10*log10(abs(GainA)),'r-');
grid on;
xlabel('Frequency (kHz)');
ylabel('Gain (dB)');
title('Acoustic efficiency (sound intensity/electrical power)');
axis([0 fMax 0 40]);