% Zrs 


clear all;
close all;

%Abrimos los ficheros de datos
flow0=load('Flow.txt');
pres0=load('Pressure.txt');

% Definimos las constantes del programa
n = size(flow0,1); % Longitud de los datos
fs = 128; % frecuencia de muestreo
t = [0:1/fs:(size(flow0,1)-1)*(1/fs)]; % Tiempo de adquisicion de datos
fp = 2; % frecuencia del filtro

% Dibujamos las seņales
figure(1),
subplot(2,1,1), plot(t,flow0);
title('flow & pressure recordings');
ylabel('Flow (l/s)');
subplot(2,1,2), plot(t,pres0);
ylabel('Pressure (hPa)');
xlabel('Time (s)');

% Calculamos la Transformada de Fourier
fflow = fft(flow0);
fpres = fft(pres0);

% Calculamos el espectro de potencia
pflow = fflow.*conj(fflow); 
ppres = fpres.*conj(fpres);

f = fs*(0:n/2)./n; % Normalizamos las unidades de frecuencia

% Representamos los datos del espectro
figure(2)
subplot(2,1,1), semilogy(f,pflow(1:n/2+1,1),'.');
title('power spectrum');
ylabel('Flow');
subplot(2,1,2), semilogy(f,ppres(1:n/2+1,1),'.');
ylabel('Pressure');
xlabel('Frequency (Hz)');

% Restramos el valor medio de las seņales
flow1 = flow0 - mean(flow0);
pres1 = pres0 - mean(pres0);

% Filtramos con el filtro 'Moving Average' para quedarnos con las
% frecuencias bajas
np = fs/fp;  % Numero de puntos de la ventana del filtro
b = (1/np)*ones(1,np);  % Numerador de la funcion de transferencia
a = 1; % Denominador de la funcion de transferencia

% Representamos el filtro
figure(3)
freqz(b,a);

% Filtramos la seņal con la funcion 'filtfilt'
fil_flow = filtfilt(b,a,flow1);
fil_pres = filtfilt(b,a,pres1);

% Representamos la seņal filtrada
figure(4)
subplot(2,1,1), plot(t,fil_flow);
title('Breathing siganals');
ylabel('Flow (l/s)');
subplot(2,1,2), plot(t,fil_pres);
ylabel('Pressure (hPa)');
xlabel('Time (s)');

% Obtenemos la seņal oscilatoria restando la seņal filtrada de la seņal sin
% filtrar
flow1 = flow1 - fil_flow;
pres1 = pres1 - fil_pres;

% Recortamos el primer y ultimo 0.5 s
flow1 = flow1((fs/fp + 1):(length(flow0)- fs/fp));
pres1 = pres1((fs/fp + 1):(length(flow0)- fs/fp));
t = [0:1/fs:(size(flow1,1)-1)*(1/fs)]; % Redefinimos el tiempo

% Representamos la seņal oscilatoria
figure(5)
subplot(2,1,1), plot(t,flow1);
title('Oscillatory signals');
ylabel('Flow (l/s)');
subplot(2,1,2), plot(t,pres1);
ylabel('Pressure (hPa)');
xlabel('Time (s)');

flow2 = flow1 - mean(flow1);
pres2 = pres1 - mean(pres1);

% Dividimos los datos en n bloques. Ponemos los bloques en una
% matriz donde cada columna es un bloque distinto.
blocs = input('Numero de bloques? (maximo 62) ','s'); % Pedimos el numero de bloques
blocs = str2num(blocs);
grp  = fix(length(flow1)/(fs/fp)); 
grp1 = fix(grp/blocs);
flow = flow1(1:(blocs*grp1*(fs/fp)));
pres = pres1(1:(blocs*grp1*(fs/fp)));
flow = reshape(flow,length(flow)/blocs,blocs);
pres = reshape(pres,length(pres)/blocs,blocs);
t = [0:1/fs:(size(flow,1)-1)*(1/fs)]; % Redefinimos el tiempo

% Restamos el valor medio ded cada bloque
meanFlow = [];
meanPres = [];
mflow = mean(flow);
mpres = mean(pres);
for i=1:length(flow)
    meanFlow = [meanFlow ; mflow];
    meanPres = [meanPres ; mpres];
end
flow = flow - meanFlow;
pres = pres - meanPres;

% Aplicamos una ventana de Hamming
% La funcion hamming(n) crea una ventana de Hamming de n puntos, y la
% funcion 'meshgrid' crea una matriz donde cada columna es la misma
% ventana repetida, para poder aplicarla a la matriz de datos.
H = input('Aplicamos ventana de Hamming? (s n) ','s');

switch H
    case 's'
        [H1 H2] = meshgrid(ones(blocs,1),hamming(size(flow,1)));
        
        flow = flow.*H2;
        pres = pres.*H2;
        
        flow2 = flow2.*hamming(length(flow2));
        pres2 = pres2.*hamming(length(pres2));

        % Dibujamos la ventana de Hamming
        figure(6)
        plot(H2(:,1));
        title('Hamming window');
        ylabel('weight');
        xlabel('index');
end

% Dibujamos los datos (solo el primer bloque)
figure(7),
subplot(2,1,1), plot(t,flow(:,1));
title('First block');
ylabel('Flow (l/s)');
subplot(2,1,2), plot(t,pres(:,1));
ylabel('Pressure (hPa)');
xlabel('Time (s)');

% Calculamos la Transformada de Fourier
n = size(flow,1); % Longitud de los datos
fflow1 = fft(flow);
fpres1 = fft(pres);

% Calculamos el expectro de energia
pflow1 = fflow1.*conj(fflow1); 
ppres1 = fpres1.*conj(fpres1);

f = fs*(0:n/2)./n; % Normalizamos la frecuencia

% Representamos los datos del espectro
figure(8)
subplot(2,1,1), semilogy(f,pflow1(1:n/2+1,1),'.');
title('First block power spectra');
ylabel('Flow');
subplot(2,1,2), semilogy(f,ppres1(1:n/2+1,1),'.');
xlabel('Frequency (Hz)');
ylabel('Pressure');

%  Calculamos la impedancia (z) i la representamos solo en los puntos que
%  nos interesan (multiples de 2 Hz).
z = fft(pres)./fft(flow);
val = find((f/2 - floor(f/2) == 0) & (f < 33));
val(1) = [];

% Calculamos las frequencias de excitacion
zs = fft(pres2)./fft(flow2);
f2 = fs*(0:length(flow2)/2)./length(flow2);
val2 = find((f2/2 - floor(f2/2) == 0) & (f2 < 33));
val2(1) = [];

% Calculamos la impedancia mediante el 'Cross Spectrum'
cros_flow = fft(flow);
cros_pres = fft(pres);

pp = cros_pres .* conj(cros_pres); 
pf = cros_flow .* conj(cros_pres);
ff = cros_flow .* conj(cros_flow);

% Calculamos los cros-espectros
Gpp = mean(pp,2);
Gpf = mean(pf,2);
Gff = mean(ff,2);

cros_z = Gpp ./ Gpf; % Calculamos la impedancia
gamma2 = (abs(Gpf).^2) ./ (Gpp .* Gff); % Calculamos Gamma^2

% Calculamos los errores
% warning off MATLAB:divideByZero
cros_zstdR = real(cros_z).*sqrt(((1-gamma2)./(gamma2*2*blocs))).*sqrt(1+(imag(cros_z)./real(cros_z)).^2);
cros_zstdI = imag(cros_z).*sqrt(((1-gamma2)./(gamma2*2*blocs))).*sqrt(1+(real(cros_z)./imag(cros_z)).^2);
val = find((f/2 - floor(f/2) == 0) & (f < 33));
val(1) = [];

% Representamos parte real y parte imaginaria de la impedancia
figure(9)
subplot(2,1,1), plot(f2(val2),real(zs(val2)),'-o'); 
hold on;
errorbar(f(val),real(cros_z(val)),cros_zstdI(val),'-or');
title('Zrs');
xlabel('Frequency (Hz)');
ylabel('Resistance(hPa·s/l)');
subplot(2,1,2), plot(f2(val2),imag(zs(val2)),'-o');
hold on;
errorbar(f(val),imag(cros_z(val)),cros_zstdI(val),'-or');
xlabel('Frequency (Hz)');
ylabel('Reactance (hPa·s/l)');

% Representamos gamma2
figure(10)
plot(f(val),gamma2(val),'.');
xlabel('Frequency (Hz)');
ylabel('Coherence');