load('C:\Users\d035190\Documents\MATLAB\namefo\nemefoPreparazioneDati\variabiliNemefo20101004-20140505DaTePrWiCorrette.mat', 'Rain')
X=Rain;
% outlieres detection from raw sensored output 
rainprocess=X;
rdif=find(abs(diff(rainprocess))>20);
figure
plot(rainprocess)
hold on
plot(rdif,rainprocess(rdif),'r^','linewidth',3)
rainprocess(rdif)=NaN;
%rainprocess(rdiff)=[];
figure 
plot(Rain)
hold on
plot(rainprocess)

% interpolating the missng terms by cubic interpolation 

X=rainprocess;
X(isnan(X)) = interp1(find(~isnan(X)), X(~isnan(X)), find(isnan(X)),'cubic');
RainInt=X;

figure
plot(RainInt);
hold on
plot(rainprocess);

% filtering ( applied moving average and savitzky golay but savitzkey golay
% is selected but aftere filtering we removed the negative mm values 
% study rain fall in mm and if -ve values 4exist or not as i am deleting it 

% filteredrainsmooth=smooth(RainFilt)
 abs(fft(RainInt))
 plot(abs(fft(RainInt)))
figure
plot(abs(fft(RainInt)))
 X_mags=abs(fft(RainInt));

% to display it in normalized range 
num_bins=length(X_mags)
figure
 plot([0:1/(num_bins/2 -1):1], X_mags(1:num_bins/2))

X_mags=abs(fft(X1));
 
 
 
[b,a]=butter(2,0.1);
[d,c]=butter(2,0.2);
[f,e]=butter(2,0.00001,'high') 
 X=H;
X(isnan(X)) = interp1(find(~isnan(X)), X(~isnan(X)), find(isnan(X)),'cubic');
X=filtfilt(b,a,X);
 

[b,a]=butter(2,0.1)
X=RainInt;
X1=filtfilt(b,a,X);

figure 
 plot(RainInt);
 hold on
 plot(X)
 
% v v bad effect 
% [f,e]=butter(2,0.00001,'high') 
 X=RainInt;
 X1=filtfilt(f,e,X);


 [f,e]=butter(2,0.5,'low') 
 X=RainInt;
 X2=filtfilt(f,e,X);
 hold on
 plot(X2)
 
 
 [g,h]=butter(2,0.3,'low') 
 X=RainInt;
 X3=filtfilt(g,h,X);
 hold on
 plot(X3)
 
 [i,j]=butter(2,0.35,'low') 
 X=RainInt;
 X4=filtfilt(i,j,X);
 hold on
 plot(X4)
 
  save lowpassfilteredrainwithorder2freq35.mat

 
 [i,j]=butter(3,0.4,'low') 
 X=RainInt;
 X4=filtfilt(i,j,X);
 hold on
 plot(X4)
 
 
 
 
 
 save lowpassfilteredrainwithorder3freq4.mat
 
  
  [i,j]=butter(3,0.01,'low') 
 X=RainInt;
 X4=filtfilt(i,j,X);
 hold on
 plot(X4)
 
 
 
%Z = smooth(Y) uses the moving average method with span 5 and
% filteredrainsmooth=smooth(RainFilt);
% X_smooth=abs(fft(filteredrainsmooth));
% num_bins2=length(X_smooth);
% figure
% plot([0:1/(num_bins2/2 -1):1], X_smooth(1:num_bins2/2));

%Z = smooth(Y,'sgolay') uses the Savitzky-Golay method with DEGREE=2,SPAN = 5, X = 1:length(Y).
    
filteredraingolay=smooth(RainInt,'sgolay');
X_sgolay=abs(fft(filteredraingolay));
num_bins3=length(X_sgolay);
figure
plot([0:1/(num_bins3/2 -1):1], X_sgolay(1:num_bins3/2));
plot(filteredraingolay)

% moving average 7 span 

% filteredrainsmooth7=smooth(RainFilt,7);
% butter worth filter design on interpolated signal 

% X_rain=abs(fft(RainInt));
% num_binsbutter=length(X_rain);
% figure
% plot([0:1/(num_binsbutter/2 -1):1], X_rain(1:num_binsbutter/2));
% [b a]= butter(2,0.1,'low');
% H=freqz(b,a, floor(num_binsbutter/2));
% hold on
% plot([0:1/(num_binsbutter/2 -1):1], abs(H),'r');
% X_filtered=filter(b,a,RainInt)
% figure
% plot(RainInt)
% hold on
% plot(X_filtered)

plot(filteredrainsmooth)
hold on
plot(filteredraingolay)
hold on
plot(filteredrainsmooth7)


% finding the negatives ( as by applying the svatzky golay filter the
% values goes in negative so qwe should delet those....

id=find(filteredraingolay<0)
figure , plot(filteredraingolay(id))
filteredraingolay(id)=0;

% %outleir 
% rainprocess=RainFilt;
% rdif=find(abs(diff(rainprocess))>20);
% figure
% plot(rainprocess)
% hold on
% plot(rdif,rainprocess(rdif),'r^','linewidth',3)

%RAIN PREPROCESSING
minval=min(filteredraingolay);
maxval=max(filteredraingolay);
norm_rain = (filteredraingolay - minval) / ( maxval - minval );

figure
plot(norm_rain)
hold on 
save('rainfiltereddata06-08-16.mat','Rain','RainInt','filteredraingolay')


% results are not thqat much good so we are going to do some more new
% experiments i.e first we will apply savitzkey golay filter for noise
% removal and than we will apply median foiloter for further smoothining
% and then we predict the rain. 

load('C:\Users\d035190\Documents\MATLAB\namefo\nemefoPreparazioneDati\variabiliNemefo20101004-20140505DaTePrWiCorrette.mat', 'Rain')
X=Rain;
% outlieres detection from raw sensored output 
rainprocess=X;
rdif=find(abs(diff(rainprocess))>20);
figure
plot(rainprocess)
hold on
plot(rdif,rainprocess(rdif),'r^','linewidth',3)
rainprocess(rdif)=NaN;
%rainprocess(rdiff)=[];
figure 
plot(Rain)
hold on
plot(rainprocess)

% interpolating the missng terms by cubic interpolation 

X=rainprocess;
X(isnan(X)) = interp1(find(~isnan(X)), X(~isnan(X)), find(isnan(X)),'cubic');
RainInt=X;

figure
plot(RainInt);
hold on
plot(rainprocess);


filteredrainsgolay = sgolayfilt(RainInt,2,5);
figure
plot(RainInt);
hold on
plot(filteredrainsgolay);

filteredrainmed=medfilt1(filteredrainsgolay);
plot(filteredrainmed);

% filtering rain with moving average filtyer 
load('C:\Users\d035190\Documents\MATLAB\namefo\nemefoPreparazioneDati\variabiliNemefo20101004-20140505DaTePrWiCorrette.mat', 'Rain')
X=Rain;
% outlieres detection from raw sensored output 
rainprocess=X;
rdif=find(abs(diff(rainprocess))>20);
figure
plot(rainprocess)
hold on
plot(rdif,rainprocess(rdif),'r^','linewidth',3)
rainprocess(rdif)=NaN;
%rainprocess(rdiff)=[];
figure 
plot(Rain)
hold on
plot(rainprocess)

% interpolating the missng terms by cubic interpolation 

X=rainprocess;
X(isnan(X)) = interp1(find(~isnan(X)), X(~isnan(X)), find(isnan(X)),'cubic');
RainInt=X;

figure
plot(RainInt);
hold on
plot(rainprocess);


% moving average filter 

windowSize = 3;

b = (1/windowSize)*ones(1,windowSize);
a=1;
y = filter(b,a,RainInt);
figure
plot(RainInt);
hold on
plot(y);

windowSize2 = 3;
b = (1/windowSize2)*ones(1,windowSize2);
a=1;
y2 = filter(b,a,RainInt);
hold on 
plot(y2);
mvfrain=y;

save('movingavgrainfilterdata.mat','Rain','RainInt','mvfrain');

% calculating error of filtered data 

ey=mean(RainInt-y);
esavitzkey=mean(RainInt-filteredrainsgolay);
emedian=mean(RainInt-filteredrainmed);

err = immse(RainInt,y);
err2 = immse(RainInt,filteredrainsgolay);
err3 = immse(RainInt,filteredrainmed);
err4= immse(RainInt, X4);