ECG-Kit 1.0
(3,531 bytes)
%% (Internal) Calculates the wavelet transform
%
% Prototype:
% ----------
%
% wt = qs_wt(x, scales, fs, q_filters)
%
% Description:
% ------------
% Calculates the Wavelet transform on signal "x", for the scales defined in
% "scales" at sampling frequency "fs". As a result the "wt" signal is
% returned.
% Optionally the "q_filters" object can be defined externally to perform
% the design of the filters only one time.
%
% Arguments :
% ----------
% x: k Signal/s of length n. A matrix of n by k
% scales: A vector indicating which scales to calculate. Default 1:4
% fs: Sampling frequency. Default 250 Hz
% q_filters: Optional filter objects designed with
% "qs_filter_design".
%
% Examples :
% ----------
%
% % Wavelet transform filter bank design
% scales = 3:5;
% sampling_rate = 512;
% CantScales = length(scales);
% MaxScales = max(scales);
% scale_idx = nan(MaxScales,1);
%
% % Design the filters externally only once.
% filters_cache_filename = ['wt_filters_' num2str(MaxScales) ' scales_' num2str(sampling_rate) ' Hz.mat' ];
% if( exist(filters_cache_filename, 'file') )
% load( filters_cache_filename );
% else
% q_filters = qs_filter_design(MaxScales, sampling_rate);
% end
%
% % Generate a mapping variable to identify the wavelet scales in case not
% % all scales are needed
% for ii = 1:CantScales
% scale_idx(scales(ii)) = ii;
% end
%
% wtECG = qs_wt(ECG, scales, sampling_rate, q_filters);
%
% % All signals at scale 4 are accesed like this.
% wtECG(:,:,scale_idx(4))
%
%
% Author: Mariano Llamedo Soria (llamedom at {electron.frba.utn.edu.ar || unizar.es})
% Birthdate: 17/2/11
% Last update: 20/2/13
function wt = qs_wt(x, scales, fs, q_filters)
if( nargin < 3 || isempty(fs) )
fs = 250;
warning('Assuming fs=250Hz.');
end
if( nargin < 2 || isempty(scales) )
scales = 1:4;
warning('Assuming 4 scales of wavelet decomposition.');
end
cant_scales = length(scales);
max_scale = max(scales);
if( nargin < 4 || isempty(q_filters) )
filters_cache_filename = ['wt_filters_' num2str(max_scale) ' scales_' num2str(fs) ' Hz.mat' ];
if( exist(filters_cache_filename, 'file') )
load( filters_cache_filename );
else
q_filters = qs_filter_design(scales, fs);
if( usejava('desktop') )
disp('Please check the filters high-pass frequency response. Press Run (F5) to save these filters for future use.')
fvtool(q_filters, 'Fs', fs)
keyboard
end
this_path = mfilename('fullpath');
this_path = fileparts(this_path);
save([this_path filesep filters_cache_filename], 'q_filters');
end
end
[CantSamp, CantSig] = size(x);
wt = zeros(CantSamp, CantSig, cant_scales);
% scale_idx = nan(max_scale,1);
for ii = 1:cant_scales
%indice para saber que escala corresponde a cada columna de MyWavelets
% scale_idx(scales(ii)) = ii;
aux_delay = grpdelay(q_filters(ii),1);
cant_casc = length(aux_delay);
aux_delay = round( sum(aux_delay) );
aux_filter_length = 2*aux_delay+cant_casc;
aux = filter(q_filters(ii), x);
wt(:,:,ii) = [ zeros(aux_delay,CantSig); aux(aux_filter_length:end, :); zeros(aux_delay+1,CantSig) ];
end