% Filterung der Messdaten
% Stefan Wakolbinger 18.03.2011
clear all; close all; clc;

load aData_10_50_7700_1.mat;

len = length(aData);    % Anzahl der Messungen
j = 1;                  % Laufvariable für subplot
rows = 3; columns = 3;  % Subplots pro Fenster
cutoff = 50;
for i=1:len
    x_sense = aData(i,1).scheibe_x;
    y_sense = aData(i,1).scheibe_y;
    
    %Filterparameter wählen
    n = round(aData(i,1).n_mean);
    
    if (n<=750)
        f_end = 0.1;
    elseif (n>7500 && n<=1400)
        f_end = 0.05;
    elseif n>1400 && n<2200
        f_end = 0.010;        
    else
        f_end = 0.0025;
    end
    
    %Filtern
    B=fir1(10,f_end);
    x_sense_filt=filtfilt(B,1,x_sense); %zero-phase filtering
    y_sense_filt=filtfilt(B,1,y_sense); %zero-phase filtering
    
%     h=fdesign.lowpass('Fp,Fst,Ap,Ast',0.001,0.02,1,60);
%     d=design(h,'equiripple'); %Lowpass FIR filter
%     x_sense_filt=filtfilt(d.Numerator,1,x_sense); %zero-phase filtering
%     y_sense_filt=filtfilt(d.Numerator,1,y_sense); %zero-phase filtering

    if mod(j,rows*columns)==1
        j=1; 
        fig=figure();
        set(fig,'units','normalized','outerposition',[0 0 1 1])
    end
    subplot(rows,columns,j); j=j+1;
    plot(1000*x_sense, 1000*y_sense,'Color',[1,0,0]);    
    plotRange=cutoff:aData(i,1).nrPoints-cutoff;
    hold on;
    plot(1000*x_sense_filt(plotRange), 1000*y_sense_filt(plotRange))
    xlabel('x-Achse in \mum');
    ylabel('y-Achse in \mum');
    title([num2str(n),' U/min'])
end

%% FFT X, von der letztgeplotten Messung
close all;
L = length(x_sense);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
X1 = fft(x_sense_filt,NFFT)/L;
X2 = fft(x_sense,NFFT)/L;
f = aData(1,1).SampleRate/2*linspace(0,1,NFFT/2);

% Plot single-sided amplitude spectrum.
plot(f,[2*abs(X2(1:NFFT/2)), 2*abs(X1(1:NFFT/2))]) 
title('Single-Sided Amplitude Spectrum of x(t)')
xlabel('Frequency (Hz)')
ylabel('|X(f)|')
legend('ungefiltert','gefiltert')

%% FFT Y, von der letztgeplotten Messung
close all;
L = length(y_sense);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y1 = fft(y_sense_filt,NFFT)/L;
Y2 = fft(y_sense,NFFT)/L;
f = aData(1,1).SampleRate/2*linspace(0,1,NFFT/2);

% Plot single-sided amplitude spectrum.
plot(f,[2*abs(Y2(1:NFFT/2)), 2*abs(Y1(1:NFFT/2))]) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
legend('ungefiltert','gefiltert')