Skip to content
Snippets Groups Projects

Introduce measurements part of the codebase

Merged Imported Introduce measurements part of the codebase
Sweeps into main
Merged Imported Achilles Kappis requested to merge Sweeps into main
Compare
and
4 files
+ 223
− 1
Compare changes
  • Side-by-side
  • Inline
Files
4
Signal Processing/Measurements/invSweep.m 0 → 100755
+ 213
− 0
%% Function to calculate the coherence function between two vectors
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 21/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate the transfer function from an Exponential Sine
% Sweep measurement.
% --------------------------------------------------
% Input arguments
%
% rec [numeric]: This can be a vector holding a single channel recording or
% a matrix in which case each column will represent a
% channel to be deconvolved.
%
% ref [numeric]: This can be a vector holding a single reference signal or
% a matrix holding a reference signal for each channel to be
% deconvolved. In case this is a matrix the number of
% columns must match the number of columns in "rec".
%
% equalise [logical] (Optional) : Whether to equalise the inverse signal or
% not to compensate for the "pink"
% spectrum. [Default: false].
%
% regFreq [numeric] (Optional): The frequency extrema for the
% regularisation process. This must be a
% vector with two real values specifying the
% frequency interval boundaries. The
% frequencies outside the specified interval
% will be highly regularised, while those
% inside the interval will be regularised
% with lower values. The regularisation
% values are defined by the "regVals"
% argument. If the sampling frequency is NOT
% specified ("fs" argument) the numbers
% provided here will correspond to frequency
% bin indices instead of frequencies.
% [Default: 0, length(reference) or 0, fs].
%
% regVals [numeric] (Optional): The regularisation values to be used. This
% must be a real vector with two elements
% specifying the in- and out-of-interval
% regularisation values to be applied.
% [Default: 0, 0].
%
% fs [numeric] (Optional): The sampling frequency. It is used mainly to
% calculate the frequency bins corresponding to
% frequencies to be regularised. [Default: NaN]
%
% causIrLen [numeric] (Optional): The lenght of the returned causal impulse
% response. This must be a a real positive
% scalar less than or equal to the length
% of the recorded signal(s).
%
% --------------------------------------------------
% Output arguments
%
% invFiltFd [numeric]: The frequency response of the inverse filter.
%
% fr [numeric]: The calculated frequency response(s).
%
% ir [numeric]: The complete (both causal and anticausal components)
% impulse after deconvolution.
%
% causIr [numeric]: The causal part of the impulse response.
%
% invFiltTd [numeric]: The impulse response of the inverse filter.
%
% --------------------------------------------------
% Notes
%
% For more information refer to "Advancements in impulse response
% measurements by sine sweeps" by Angelo Farina.
%
% --------------------------------------------------
function [invFiltFd, fr, ir, causIr, invFiltTd] = invSweep(rec, ref, equalise, regFreq, regVals, fs, causIrLen)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(2, 7); % Check for number of arguments
nargoutchk(0, 5); % Check for number of arguments
% ====================================================
% Validate input arguments
% ====================================================
% Validate mandatory arguments
validateattributes(rec, {'numeric'}, {'2d', 'nonempty', 'real', 'finite', 'nonnan'}, mfilename, "Recorded signal", 1);
if isvector(rec)
validateattributes(ref, {'numeric'}, {'vector', 'nonempty', 'real', 'finite', 'nonnan'}, mfilename, "Reference signal", 2);
rec = rec(:);
ref = ref(:);
else
validateattributes(ref, {'numeric'}, {'2d', 'nonempty', 'real', 'finite', 'nonnan'}, mfilename, "Reference signal", 2);
if ~isvector(ref) && size(ref, 2) ~= size(rec, 2)
error("The number of columns in the recorder and reference signals must be the same");
elseif isvector(ref)
ref = ref(:);
end
end
% Validate optional arguments
if nargin > 2 && ~isempty(equalise)
validateattributes(equalise, {'logical'}, {'scalar', 'nonnan', 'finite'}, mfilename, "Equalise for pink spectrum", 3);
else
equalise = false;
end
if nargin > 3 && ~isempty(regFreq)
if isempty(fs) || isnan(fs) || ~isfinite(fs)
validateattributes(regFreq, {'numeric'}, {'vector', 'real', 'finite', 'nonnan', 'nonnegative', 'increasing', 'integer', 'numel', 2}, mfilename, "Regularisation interval extremum frequencies", 4);
else
validateattributes(regFreq, {'numeric'}, {'vector', 'real', 'finite', 'nonnan', 'nonnegative', 'increasing', 'numel', 2}, mfilename, "Regularisation interval extremum frequencies", 4);
end
else
regFreq = NaN;
end
if nargin > 4 && ~isempty(regVals)
validateattributes(regVals, {'numeric'}, {'vector', 'real', 'finite', 'nonnan', 'nonnegative', 'numel', 2}, mfilename, "Regularisation values", 5);
else
regVals = zeros(2, 1);
end
if nargin > 5 && ~isempty(fs)
validateattributes(fs, {'numeric'}, {'scalar', 'real', 'finite', 'nonnan', 'positive'}, mfilename, "Sampling frequency", 6);
else
fs = NaN;
end
if nargin > 6 && ~isempty(causIrLen)
validateattributes(causIrLen, "numeric", {'scalar', 'real', 'finite', 'nonnan', 'nonempty', 'positive'}, mfilename, "The length of the causal impulse response to be returned", 7);
if isvector(rec) && causIrLen > numel(rec)
error("Returned causal impulse response cannot be larger than the provided recorded signal");
elseif ~isvector(rec) && causIrLen > size(rec, 1)
error("Returned causal impulse response cannot be larger than the provided recorded signals");
end
else
causIrLen = size(rec, 1);
end
% ====================================================
% Calculate inverse filters
% ====================================================
invFiltFd = fft(ref, size(rec, 1) + size(ref, 1) - 1);
% Check for cases
if ~equalise
invFiltFd = conj(invFiltFd);
elseif equalise && isnan(regFreq)
% The formula below is the so called "Kirkeby Inverse Filter"
% presented in Farina's paper for frequency domain inversion.
% Here it is used with a constant regularisation value.
invFiltFd = conj(invFiltFd)./(conj(invFiltFd) .* invFiltFd + 1e3 * eps);
else
% Get the bins for regularisation
if isnan(fs)
freqBins = 1:ceil(numel(invFiltFd)/2);
else
% Calculate the frequency bins for the used FFT length
freqBins = fs * (0:1/numel(invFiltFd):1/2 - 1/numel(invFiltFd));
end
% Get regularisation values for each bin
lRegIdx = find(freqBins < regFreq(1), 1, 'last'); % Low frequency bin index
hRegIdx = find(freqBins > regFreq(2), 1, 'first'); % High frequency bin index
% Make sure the frequency indices are valid
if isempty(lRegIdx) || isempty(hRegIdx)
if ~isnan(fs)
error("The regularisation interval frequencies provided are not valid.")
else
error("The regularisation interval bin indices provided are not valid.")
end
end
% Create regularisation values vector
regFacs = ones(size(freqBins)) * regVals(2); % High regularisation values (out-of-band)
regFacs(lRegIdx:hRegIdx) = regVals(1); % Low regularisation values (in-band)
regFacs = regFacs(:); regFacs = [regFacs; flip(regFacs(2:end))];
% Calculate the frequency-dependent regularised inverse filter
invFiltFd = conj(invFiltFd)./(conj(invFiltFd) .* invFiltFd + regFacs);
end
% ====================================================
% Convolve with inverse filter to get frequency response
% ====================================================
if nargout > 1
recFFT = fft(rec, numel(invFiltFd)); % Get the frequency response of the recorded signal
fr = recFFT .* invFiltFd; % Perform deconvolution
end
if nargout > 2
ir = ifft(recFFT .* invFiltFd); % Calculate impulse response
end
if nargout > 3
causIr = ir(1:causIrLen, :);
end
if nargout > 4
invFiltTd = ifft(invFiltFd);
end
end
\ No newline at end of file