function pulsanalyse()

figure(1);clf
x=csvread('photopl.txt',4,0);
x=x/std(x);
t=((0:length(x)-1)/40)';
plot(t,x)
xlabel('t [sek]')
ylabel('x(t)')

n=findpeaks(x);
hold on
plot(t(n),x(n),'go')
axis([0 t(end) -2 5])
xlabel('t [sek]')
ylabel('x(t)')

figure(2);clf
window=hanning(256);
noverlap=0;
nfft=1024;
fs=40;
[Pxx,f] = pwelch(x,window,noverlap,nfft,fs);
subplot(2,1,1)
plot(f,10*log10(Pxx))
grid
xlabel('f [Hz]')
ylabel('Magnitude [dB]')
axis([0 20 -40 20])

b=fir1(150,[0.58 1.27]/(40/2));
[h,w]=freqz(b);
subplot(2,1,2)
plot(w/(2*pi)*40,20*log10(abs(h)));
grid
xlabel('f [Hz]')
ylabel('Magnitude [dB]')
axis([0 20 -100 20])

y=filtfilt(b,1,x);
figure(1)
hold on
plot(t,y,'r')
n2=findpeaks(y);
plot(t(n2),y(n2),'ko')

figure(2)
[Pyy,f] = pwelch(y,window,noverlap,nfft,fs);
subplot(2,1,1)
hold on
plot(f,10*log10(Pyy),'r')
axis([0 20 -100 20])

d=diff(t(n2));
bpm=1/mean(d/60);
rrvar=std(d);

