function [q1, q2, q3, q4,q5, messages] = qspfilt5(fs,messages)
% This function obtains the quadratic splines wavelet filterbank filters from
% scale 1 to 4 as a function of the sampling frequency, in order to use filters
% with similar analog frequency behaviour for diferent sampling frecuencies.
%
% Juan Pablo Martínez Cortés
%
% Last update: Rute Almeida 24/11/11
%
% Designed for MATLAB Version R12; tested with MATLAB Version R13
%
if nargin<2
messages = new_empty_mesg;
end
messages.status=1;
% beutors h(n) y g(n)
h = 0.125*[1 3 3 1]; %high pass % paso baxo
g = 2*[1 -1]; %low pass % paso alto
h2 = kron(h,[1 0]);
g2 = kron(g,[1 0]);
h4 = kron(h,[1 zeros(1,3)]);
g4 = kron(g,[1 zeros(1,3)]);
h8 = kron(h,[1 zeros(1,7)]);
g8 = kron(g,[1 zeros(1,7)]);
g16=kron(g,[1 zeros(1,15)]); % Rute 21/05/02 %%
q1 = g; % Q1(w) = G(w);
q2 = conv(h,g2); % Q2(w) = H(w)G(2w);
q3 = conv(conv(h,h2),g4); % Q3(w) = H(w)H(2w)G(4w);
q4 = conv(conv(conv(h,h2),h4),g8); % Q4(w) = H(w)H(2w)H(4w)G(8w);
q5 = conv(conv(conv(conv(h,h2),h4),h8),g16); % Q5(w) = H(w)H(2w)H(4w)H(8w)G(16w); % Rute 21/05/02 %%
l1 = max(find(q1~=0)); %#ok<MXFND>
l2 = max(find(q2~=0)); %#ok<MXFND>
l3 = max(find(q3~=0)); %#ok<MXFND>
l4 = max(find(q4~=0)); %#ok<MXFND>
l5 = max(find(q5~=0)); %#ok<MXFND> % Rute 21/05/02 %%
q1 = q1(1:l1)';
q2 = q2(1:l2)';
q3 = q3(1:l3)';
q4 = q4(1:l4)';
q5 = q5(1:l5)'; % Rute 21/05/02 %%
% nao foi alterado para adaptar q5 a outras frequencias de amostragem ainda
% Adautazión de os filtros de frec. 250 Hz ta atras frecuenzias de muestreo
% Adaptation of filters from freq. 250 Hz to other sampling frequencies
% 500 Hz
if fs~=250,
if fs == 500,
q1 = interp([0;q1;0],2,1); q1 = q1(2:end-2);
q2 = interp([0;q2;0],2,3); q2 = q2(2:end-2);
q3 = interp([0;q3;0],2,7); q3 = q3(2:end-2);
q4 = interp([0;q4;0],2,7); q4 = q4(2:end-2);
q5 = interp([0;q5;0],2,7); q5 = q5(2:end-2);
elseif (fs == 1000),
q1 = interp([0;q1;0],4,1); q1 = q1(2:end-4);
q2 = interp([0;q2;0],4,3); q2 = q2(2:end-4);
q3 = interp([0;q3;0],4,7); q3 = q3(2:end-4);
q4 = interp([0;q4;0],4,7); q4 = q4(2:end-4);
q5 = interp([0;q5;0],4,7); q5 = q5(2:end-4);
elseif fs == 200,
%q1 = interp([0;q1;0],4,1,0.4); q1 = q1(2:end-4);
q2 = interp([0;q2;0],4,3,0.4); q2 = q2(2:end-4);
q3 = interp([0;q3;0],4,7,0.4); q3 = q3(2:end-4);
q4 = interp([0;q4;0],4,7,0.4); q4 = q4(2:end-4);
q5 = interp([0;q5;0],4,7,0.4); q5 = q5(2:end-4);
% No cal filtrar-ie pos ya ye filtrau a 100 Hz.
%q1 = q1(1:5:end);
q2 = q2(4:5:end);
q3 = q3(5:5:end);
q4 = q4(2:5:end);
q5 = q5(2:5:end);
q1 = 1.1*[5/4, -5/4];
elseif fs == 360,
% Interpolo por 6
q1 = interp([0;q1;0],6,1); q1 = q1(2:end-6);
q2 = interp([0;q2;0],6,3); q2 = q2(2:end-6);
q3 = interp([0;q3;0],6,7); q3 = q3(2:end-6);
q4 = interp([0;q4;0],6,7); q4 = q4(2:end-6);
q5 = interp([0;q5;0],6,7); q5 = q5(2:end-6);
% Diezmo por 5
q1=q1(4:5:end);
q2=q2(1:5:end);
q3=q3(5:5:end);
q4=q4(3:5:end);
q5=q5(3:5:end);
% Interpolo por 6
q1 = interp([0;q1;0],6,1); q1 = q1(2:end-6);
q2 = interp([0;q2;0],6,3); q2 = q2(2:end-6);
q3 = interp([0;q3;0],6,7); q3 = q3(2:end-6);
q4 = interp([0;q4;0],6,7); q4 = q4(2:end-6);
q5 = interp([0;q5;0],6,7); q5 = q5(2:end-6);
%Diezmo por 5
q1 = q1(2:5:end);
q2 = q2(5:5:end);
q3 = q3(4:5:end);
q4 = q4(4:5:end);
q5 = q5(4:5:end);
else
messages.errors=[messages.errors {'There are no wavelets designed for this sampling frequency.'}];
% warning(char(messages.errors(end)))
messages.errors_desc=[messages.errors_desc 'Default wavelet filters not defined for thia samplig frequency. Use filter design instead.'];
messages.status=0;
end
end