function [v,tv] = reedMIC(wavfile,csvfile,s,U,vmax,varargin)
% [v,tv] = reed(wavfile,csvfile,s,U,vmax) gibt den Geschwindigkeitsvektor v und
% den dazugehörigen Zeitvektor tv zurück. wavfile ist der Pfad zur
% .wav-Datei; csvfile ist der Pfad zur Ausgabedatei; s ist der
% Schwellenwert ab dem ein Peak erkannt werden soll; U ist der
% Radumfang; vmax ist die theoretisch maximale erkennbare Geschwindigkeit
% Es findet keine Filterung der Daten statt.
%
% [v,tv] = reed(...,dofilter) ist dofilter > 0 wird ein
% Bandpass mit den Grenzfrequenzen 3.2kHz und 3.5kHz angewendet.
%
% [v,tv] = reed(...,F3dB1,F3dB2) die Variablen F3dB1 und
% F3dB2 sind die untere bzw. die obere Grenzfrequenz für den Bandpass.
[data,fs]=wavread(wavfile);
dataUF = data;
if length(varargin)>0 && varargin{1}>0
F3dB1 = 3.2e3;
F3dB2 = 3.5e3;
if length(varargin)==2
F3dB1 = varargin{1};
F3dB2 = varargin{2};
end
D = fdesign.bandpass('N,F3dB1,F3dB2',2,F3dB1,F3dB2,fs);
Hd= design(D);
[B,A] = sos2tf(Hd.sosMatrix);
data = filtfilt(B,A,data).*Hd.ScaleValues(1);
end
eventlength=U*fs/vmax*3.6;
ww = eventlength/2;
if mod(ww,2) ~= 0
ww = ww -1;
end
diffdata=diff(data);
t=[1:length(diffdata)]./fs;
tevtup=t(diffdata >s);
tevtdown=t(diffdata <-s);
tevtup([0 diff(tevtup)<ww/fs]>0)=[];
tevtdown([0 diff(tevtdown)<ww/fs]>0)=[];
figure();hold on;
plot(t,diffdata,'b-')
plot(tevtup,0,'r*')
plot(tevtdown,0,'cs')
figure();hold on;
vup=[U./diff(tevtup)*3.6]';
vdown=[U./diff(tevtdown)*3.6]';
tvup=zeros(length(vup),1);
tvdown=zeros(length(vdown),1);
for i=1:length(tevtup)-1
tvup(i) = mean(tevtup(i:i+1));
end
for i=1:length(tevtdown)-1
tvdown(i) = mean(tevtdown(i:i+1));
end
plot(tvup,vup,'r-');
plot(tvdown,vdown,'c-');
if length(tvup)==length(tvdown)
tv = mean([tvup,tvdown],2);
v = mean([vup,vdown],2);
else
tv=unique([tvup;tvdown]);
v=mean([interp1(tvup,vup,tv),interp1(tvdown,vdown,tv)],2);
end
plot(tv,v,'b-');
xlabel('Zeit [s]');ylabel('Geschwindigkeit [km/h]');
export=[tv,v];
if ~isempty(csvfile)
save(csvfile,'export','-ascii','-double','-tabs');
end