function tilpass(dataFolder,saveFile)
% Utfører minste-kvadrat-tilpassing til ein sinusfunksjon på to lydsignal 
% for å finne faseforskjellen mellom dei (modulo 2\pi).
% Tar to input: dataFolder og saveFolder.
% dataFolder: en string med mappen som inneholder dataene. 
% saveFile:   en string med filen som inneholder dataene. 
%
% Eksempel på bruk: tilpass('Data1','Data1\Res\res.mat')
%
% Skrevet av Pål Ellingsen, v 16.09.2011 delvis basert på
% sinusregresjon, skrevet av Iver Bakken Sperstad, som igjen var bassert på
% tilpassingsrutine for Cavendish-lab i Mekanisk fysikk av Frantz Stabo-Eeg
% og Eirik Glimsdal (2006-2008).

files = dir([dataFolder,'\*.mat']) %finner alle filene i mappen med ending .mat

for k=1:length(files) % går gjennom alle filene og tilpasser på hver av de

    % Les data fra fil (antek ei kolonne for kvart signal)
    
    files(k).name
    filename = [dataFolder,'\',files(k).name]
    if isempty(who('-file', filename ,'-var','lydlab')) % tester at det er rett fil
        % ikke rett
        if input(['Filen ', filename,...
		' inneholder ikke måledat. Vil du fortsette til neste fil? y/n'])=='y';
            continue
        else
            return
        end
    end
  
    load(filename); % laster inn alle variablene 

    % Initialiserer tilpassingsparametrar
    % (antek signal på forma s(t) = S + A*sin(2*pi*f*t + phi)

    S = [0 0];      % middelverdi for signalet
    A = [0 0];      % amplitude
    f = [0 0];      % frekvens
    phi = [0 0];    % fase

    % Går gjennom eitt av signala om gongen
    for signal = 1:numChan

        s = s1(:,signal); % plukkar ut rett signal
        len = length(s);
        t = (1:len)';       
	% tidsvektor med heiltal; tek ikkje omsyn til samplingsfrekvens før på slutten

        % Initiell tilpassingsrutine, steg 1  
        minS = min(s);
        maxS = max(s);
        meanS = 0.5*(minS + maxS);          
	% finn omtrentleg middelverdi for signalet (dvs. gitt at talet på periodar er heiltal)

        prevAboveMean = ( s(1) > meanS );   
	% finn ut om signalet startar over (= 1) eller under (= 0) middelverdien
        crossMean = [];                     
	% kryssingspunkt, dvs. når signalet kryssar middelverdien

        % Initiell tilpassingsrutine, steg 2: Gå gjennom signal og lokaliser kryssingspunkta
        for i=2:len
            aboveMean = ( s(i) > meanS );       % er signalet over eller under middelverdi nå?
            if ( aboveMean ~= prevAboveMean )   
		% dersom signal går frå å vere over(under) til å vere under(over) middelverdi, ...
                crossMean = [crossMean i];      % ...lagre som kryssingspunkt
                prevAboveMean = aboveMean;
            end
        end

        % Initiell tilpassingsrutine, steg 3: 
	% Finn ut kor mange periodar signalet består av (med korrigering av støy i signalet)

        halfPeriod = crossMean(2:end) - crossMean(1:end-1);         
	% lengda på halvperiodar, dvs. avstand mellom kryssingspunkt 

        numHalfPeriods = length(halfPeriod);

        badHalfPeriods = find( halfPeriod < max(halfPeriod)/2 );    
	% dersom ei halvperiode er uvanleg kort skuldast det 
	% truleg berre støy rundt kryssingspunkt, ...
	
        numHalfPeriods = numHalfPeriods - length(badHalfPeriods);   
	% ...så ikkje ta omsyn til desse i berekning av periodetal

        meanPeriod = 2*sum(halfPeriod)/numHalfPeriods;        % estimerer lengda på ein periode

        if (numHalfPeriods < 3)
            disp('Feil: Klarer ikkje å finne ein heil periode i signalet.');
            return;
        end

        % Set startverdiar for tilpassingsrutine basert på initiell tilpassingsrutine
        S0 = meanS;
        A0 = maxS - meanS;           
        f0 = 1/meanPeriod;

        phi0 = -2*pi*f0*crossMean(1);                           
	% baserer fase-estimat på fyrste registrerte kryssingspunkt
        if ( s( crossMean(1) + round(0.25*meanPeriod) ) < S0 )  
	    % dersom signalet er på veg nedover (og ikkje oppover) ved fyrste kryssingspunkt...
            phi0 = phi0 - pi;                                   
	    % ...så må fasen korrigerast tilsvarande ein halv periode
        end

        x0=[S0,A0,2*pi*f0,phi0];     % startpunkt for tilpassingsiterasjonar
        
        modelFunc = @(x,t) x(1) + x(2).*sin(x(3).*t +x(4)); % lager modellfunksjon
    
        sumSquareDeviation=@(x,s,t) sum((s-modelFunc(x,t)).^2 ); % regner ut kvadratavviket
        % Kurvetilpassing ved optimalisering av tilpassingsparametrare
        [xFit fval] = fminsearch(@(x) sumSquareDeviation(x,s,t), x0);

        % Lagrar tilpassingsparametrar
        S(signal) = xFit(1);
        A(signal) = xFit(2);
        f(signal) = xFit(3)/(2*pi)*sampleFreq;  % tek omsyn til samplefrekvens her
        phi(signal) = xFit(4);     
	% NB: tilpassingsparameteren er ikkje nødvendigvis avgrensa til eitt 2\pi-intervall

        % Bereknar kurve basert på tilpassingsparametrar
        sFit=modelFunc(xFit,t);

        % Plottar signal (tidsakse skalert med samplingsfrekvens, 
        % skalerer signal slik at middelverdi blir null for båe)
        figure(1)
        if (signal == 1)
            plot(t/sampleFreq,s-S(1),'b.');
            hold on
            plot(t/sampleFreq,sFit-S(1),'r--')
        elseif (signal == 2)
            plot(t/sampleFreq,s-S(2),'Color',[0 0.5 0],'Marker','.','LineStyle','None');
            plot(t/sampleFreq,sFit-S(2),'m--')
            hold off
            xlabel('t [s]')
            ylabel('Lufttrykk (skalert og med vilkårlege einingar)')
        end

    end

    legend('Signal 1','Kurvetilpassing til signal 1','Signal 2','Kurvetilpassing til signal 2')
    
    phi=mod(phi,2*pi); % Forandrer phi slik at den blir mellom [0,2\pi)
    
    % Skriv resultat til skjerm
    disp(' ')
    disp('*** Tilpassingsparametrar for hhv. signal 1 og 2 ***')
    disp(' ')
    disp(['Amplitude: ' num2str(A)])
    disp(['Frekvens:  ' num2str(f)])
    disp(['Fase:      ' num2str(phi)])       
    
    figure(1)
    
    serBraUt = input('Er tilpassnigen bra? y/n [y]','s');
    if isempty(serBraUt) || serBraUt == 'y'
    else
        return
    end
    
    if isempty(whos('mic1')) % tester om mic1 finnes
        resultater(k,:) = [A,f,phi];
        formatResultater = 'mic1, mic2, A1 ... An, f1... fn, phi1... phin'; 
	% viser hvordan dataene er lagret
    else
        resultater(k,:) = [mic1,mic2, A,f,phi];
        formatResultater = 'mic1, mic2, A1 ... An, f1... fn, phi1... phin'; 
	% viser hvordan dataene er lagret
    end

end

save(saveFile,'resultater','formatResultater') % lagerer dataene.

end
