//datacrunchpitotr2.sce R2
//www.basicairdata.eu JLJ 2013
//Sample data handling of pitot sensor data
//Data file "hcla12x1hz 8 deg C looking for drift" must be placed in the current directory
clear all;
xdel(winsid())
//y=read("hcla12x1hz 8 deg C looking for drift",-1,1);
y=read("hcla12x1hz 20 deg C looking for drift",-1,1);
//y=read("hcla12x10Hz",-1,1);

samplerate=1;//Hz must match with file data

meany=mean(y);
sdevy=stdev(y);
[samplesize zz]=size(y)
samplet=linspace(1,samplesize,samplesize)
hand1=figure(1,"BackgroundColor",[1,1,1]);
plot(samplet,y);
title("Sensor Output plotting","fontsize",4)
xlabel("Sample number","fontsize",4)
ylabel("Decimal count, sensor output","fontsize",4)
hand2=figure(2,"BackgroundColor",[1,1,1]);
histplot(5,y,normalization=%f)
title("Histogram of samples","fontsize",4)
xlabel("Sample value, divided in bins","fontsize",4)
ylabel("Total occurences","fontsize",4)

//Remove offset z=y-meany
z=y-meany;
hand=figure(3,"BackgroundColor",[1,1,1]);
plot(samplet,z);
title("Sensor Output plotting, mean substracted","fontsize",4)
xlabel("Sample number","fontsize",4)
ylabel("Decimal count, sensor output","fontsize",4)
meanz=mean(z);
sdevz=stdev(z);
disp(meanz);
disp(sdevz);
//Linear regression
[a,b,sig]=reglin(samplet,z')
scf(hand);
regresso=a.*samplet+b
plot(samplet,regresso,'red')

//Moving average filter, output calculation
N=50; //Number of samples to average
Bm=ones(N,1); 
zsmooth=1/N*conv(z,Bm);
plot(samplet,zsmooth(1:samplesize,1),'yellow');

//FFT plot

yf=fft(zsmooth);
f=samplerate*(0:(samplesize/2))/samplesize; 
n=size(f,'*')
hand4=figure(4,"BackgroundColor",[1,1,1]);
plot(f,abs(yf(1:n)))
title("FFT smothed measurement series","fontsize",4)
xlabel("Frequency","fontsize",4)
ylabel("FFT","fontsize",4)

