%Trolldist BETA v 7.60
%Dist & frequency response measuring program using log-swept sine
%Digital Eq+digital crossover
%Realtime IIR and FFT based FIR
%Mikael Bohman VT06
%MATLAB Version 7.1.0.246 (R14) Service Pack 3
%MATLAB Version 6.5.0 (R13)

close all
clear all
clc
Trollver=7.60;  %Ver.

T=7;            %Sveptid [s] Lämpligt val är inom området 0.5-15 s. T>2*T60 dvs sveptiden bör vara större än 2*efterklangstiden för rummet som man mäter i 
AMPdB=0;       %OBS Alltid <=0 . Nivå på utsignalen vid frekvensgångsmätning relativt den möjliga maxnivån. Här kan amplituden väljas tex så att 2.83 V generas ut från förstärkaren om INITTEST körs.

%INITIALT TEST MED SINUS, är bra för att bestämma vilken utspänning man
%andvänder vid mätningen. 2.83 V motsvarar 1 elektrisk W i 8 ohms last.
FREQ=100;      %Frekvens vid initaltest [Hz]
Time=10;        %Längden på initialtest [s]

%KALIBRERINGSSEKVENS, är bra om man vill absolutkalibrera ljudtrycksnivån.
%Kräver tillgång till SPL referens. Kör INITTEST=1 samt Kalibrera=1 och mät
%ljudtrycksnivån vid mätmirkofonen vid kalibreringstonen.
CALdB=95.2;         %SPL nivå vid kalibrering [dB], skriv in vad SPL-metern visade här!
OFFSET_dB=0;   %Offsetvärdet genererat av kalibreringssekvens, Uppdatera denna efter kalibrering!
CALFRERQ=1000;      %Kalibreringsfrekvens [Hz]

%SVEPT SINUS

f_start=20;    %Startfrevens för svepet [Hz]
f_end=20000;%   %Slutfrekvens för svepet [Hz]
%VIKTIGT om du använder mik. in på ljudkortet utan att använda en LINE-IN
%som referens måste PREDELAY >0 tex 64 eller 128 beroende på fs
PREDELAY=0;     %Predelay i sample (DEFAULT=0) Använd predaly vid mätning i närfält (antivikningsfiltret på ljudkortet förringer) och om ej referenskanalen (LINE-IN) används, tex PREDELAY=128;
Prefactor=2;    %Chirpen börjar på en frekvens som är f_start/Precator ÄNDRA INTE DENNA OM DU INTE VET VAD DET INNEBÄR
Postfactor=1.1;    %Chirpen slutar på en frekvens som är f_end/Postfactor ÄNDRA INTE DENNA OM DU INTE VET VAD DET INNEBÄR
%Postwindow=200;             %1/Postwindow fönstras i slutet

%DIST. MÄTNING VID KONSTANT LJUDTRYCKSNIVÅ

Konstant_level_start=300;    %Den frekvens som konstant ljudtrycksnivå börjar på, var FÖRSIKTIG om ni mäter diskanter, de tål inte hög elektrisk ineffekt!
Konstant_level_end=8000;     %Den frekvens som konstant ljudtrycksnivå slutar på
Konstant_level_LEVEL=[90 80];     %Den önskade ljudtrycksnivånerna vid dist mätning

%Digital EQ
Comp_high=15000; %Den övre gränsfrekvensen för kompensatorn [Hz]
Comp_low=200;    %Den undre gränsfrekvensen för kompensatorn [Hz]
%DISPLAY

RES=150;                    %Antalet punkters upplösning i frekvensled på logskala
ZOOM_ms=25;               %Antalet ms som visas vid fönstringen av IR
%Här kan ni själva definiera vid vilka frekvensen som ska visas med siffror
%i THD plotten. Använd detta för att få snygga plottar.
XTICK=[100;300;1000;3000;10000];    %Frekvens i [Hz]

%Ljudkort
%VIKTIGT skriv in de samplingshastigheter som DU VET att ditt ljudkort
%stödjer UTAN att skapa dist pga undermålig interpolering. 
SOUNDCARD=[16000 22050 44100 48000];    %Samplingsfrekvenser som ditt ljudkort stödjer utan att skapa vikningsdist vid interpolering.

%Spara till fil
PATH='I:\JAVAS\JAVAS Audio\MATLAB\LILLKUB MÄTNINGAR\';  %Sökväg att spara mätningar i
STRING=sprintf('Files will be saved under %s',PATH);

%MÄT n st övertoner
MAXharmonic=3;              %Hur många övertoner beräknas totalt (Påverkar THD) 3>=MAXharmonic

%Val av fönstringstidpunkter, om ej musinput används
WINDOW=[eps 1 300 400];   %Fönstringsval [ms], om ej musinput används

%VAL [0/1]
USE_WINDOW=0;   %Använd mus-baserad fönstring = 0 / Använd ej mus-baserad fönstring = 1 
INITTEST=0;     %Initialtest av nivån avslagen = 0 / Initialtest av nivån påslagen=1
KONTROLL=1;     %Grafisk kontroll av inveren avslagen = 0 / Grafisk kontroll av inversen (rekomenderas för nybörjare) = 1
Konstant_level=0;   %Mät ej dist vid konstant utnivå = 0 / Mät dist vid konstant utnivå = 1 
Simulate=1;         %Simulera system = 0 / Mät system = 1
Kalibrera=0;    % Mät system, kalibrering redan gjord = 0 / Kalibrera enbart = 1
DISP_WINDOW=1;  %Plotta alla fönster = 1 / Plotta inte oväsentliga fönster = 0
TICK=0;         %Plotta THD med automatisk frekvensskala =0 / Plotta THD med manuell frekvensskala=1, se XTICK! 
TESTcrossover=0; %Beräkna samt exportera digital EQ =1
INSIDEBOX=0;   %Mätning i frifält =0 / Mätning innuti lådan = 1

%Högtalarelement
ELEMENT='CA18RLY in 21 liter'; %Strängen måste fungera som filnamn också, inga / _  här!
CONDITION='';%  On axis eller offaxis
VOLTAGE='2.83V';
MIKAVSTAND=50; %Mätavstånd till mikrofonen

%-----------------------------------------------------------------------
%Här börjar programmet

STRING=sprintf('Trolldist %2.2f',Trollver);
disp(STRING)

%Beräkna minsta möjliga fs;
[APA,I]=find(SOUNDCARD>2*f_end*Postfactor);
if(isempty(I))
    fs=max(SOUNDCARD);
    STRING=sprintf('Impossible to create svept sine up to the frequency %d Hz !',f_end);
    disp(STRING);
    STRING=sprintf('f_end must be < %d Hz with the avaiable soundcard',floor(fs/2/Postfactor));
    disp(STRING);
    return
end
fs=SOUNDCARD(I(1));
STRING=sprintf('The samplingrate was choosen to %d Hz',fs);
disp(STRING);

%Beräkna amplituden från nivån
AMP=10^(AMPdB/20);
ZOOM=ceil(eps+fs*ZOOM_ms/1000);

%Antivikningsskydd
if((f_end*Postfactor )>fs/2) 
    f_end=fs/2/Postfactor;
end

%Initialtest av nivån, bra för kalibrering, bra för att testa om man får
%något ljud ut.
if(INITTEST>0)  ;
    disp('Playing test signal');
    t=0:1/fs:Time;
    SIN=0.99*AMP*sin(2*pi*FREQ*t);
    wavplay(SIN(:)*[1 1],fs,'sync');
    clear SIN t DATA
end

disp('Calculating signal');
t=0:1/fs:T-1/fs;
f0=f_start/Prefactor;t1=T;f1=f_end*Postfactor;phi=-90;
instPhi = t1/log(f1/f0)*(f0*(f1/f0).^(t/t1)-f0); %Beräkna fas för chirp
y=0.99*cos(2*pi*(instPhi+phi/360)); %Beräkna chirp

clear instPhi t1 f0 f1 phi

time_k=log( Postfactor * f_end ./ f_start *Prefactor)/T;    %Beräknar tidskostanten för inversen

Prewindow=floor(0.8*fs*log(Prefactor)/time_k); % Beräknar ett optimalt förfönster
Postwindow=ceil(0.8*fs*(T-log(f_end*Prefactor/f_start)/time_k)); % Beräknar ett optimalt efterfönster


%Beräkna fönster
n=linspace(-pi,0,Prewindow);
m=linspace(0,pi,Postwindow);
y(end-length(m)+1:end)=y(end-length(m)+1:end).*0.5.*(1+cos(m));  %Fönstra slutet
y(1:length(n))=y(1:length(n)).*0.5.*(1+cos(n));                 %Fönstra början

%Beräkna inversen i tidsdomänen
f=linspace(0,fs/2,length(y));
yinv=y(end:-1:1).*exp(-t*time_k);                           %Analytisk invers i tidsdomänen
clear t
YINV=fft(yinv,2*length(y));                             %Inversen i frekvensdomänen
Y=fft(y,2*length(y));                                   %Signalen i frekvensdomänen
[APA,I]=find(and(f>f_start,f<f_end));
P=mean(abs(Y(I).*YINV(I)).^2);                 %Beräkna medeleffekten av Y(f)*YINV(f) inom frekvensintervallet [f_start f_end]
yinv=yinv./sqrt(P);                                     %Korrigera amplituden på inversen

%Grafisk kontroll av inversen
if(KONTROLL>0)
    
    dBkontroll=20*log10(abs(Y(1:end/2)));
    dBkontroll_inv=20*log10(abs(YINV(1:end/2)./sqrt(P)));
    dBprodukt=20*log10( abs( Y(1:end/2).*YINV(1:end/2)./sqrt(P)));
    figure,semilogx( f, dBkontroll_inv, f , dBkontroll, f , dBprodukt );
    xlim([f_start/Prefactor f_end*Postfactor]);
    xlabel('Frequency [Hz]');
    ylabel('Level [dB]');
    legend('Inverted signal in frequency domain','Signal in frequency domain','Product in frequency domain');
    grid on
    hold on
    plot(f_start*[1 1],ylim,'k',f_end*[1 1],ylim,'k','linewidth',3)
    hold off
    ylim([-80 80])
    drawnow
    dBERROR=max(dBprodukt(I))-min(dBprodukt(I));
    STRING=sprintf('The error from the inverse will be within +- %2.3f dB\nIncrease T to achive a smaller error!',dBERROR/2);
    disp(STRING)
end
clear Y YINV f P dBkontroll dBkontroll_inv dBprodukt APA

if(Simulate>0)
    %Mät dist
    disp('Recording signal')
    pause(0.1);
    wavplay(AMP*y(:)*[1 1],fs,'async');REC=wavrecord(length(y),fs,2);
    disp('Processing recording');
    MAXVEC=max(abs(REC));
    
    if(max(MAXVEC)>0.9)
        if(MAXVEC(1)>0.9)
            disp('!WARNING!, input signal from microphone is overloaded')
        end
        if(MAXVEC(2)>0.9)
            disp('!WARNING!, input reference signal is overloaded')
        end
        disp('!!! PROGRAM HALTED !!!')
        return
        
    end
    
    
    
    if(INSIDEBOX==0)
        DATA=REC(:,1);
    elseif(INSIDEBOX==1)
        DATA=300*diff(REC(:,1),2); %Mätning med mikrofon i lådan
    end
    REF=REC(:,2);
    DATA(end-length(m)+1:end)=DATA(end-length(m)+1:end).*0.5.*(1+cos(m)');  %Fönstra det inspelade slutet
    DATA(1:length(n))=DATA(1:length(n)).*0.5.*(1+cos(n)');                    %Fönstra den inspelade början
    REF(end-length(m)+1:end)=REF(end-length(m)+1:end).*0.5.*(1+cos(m)');  %Fönstra det inspelade slutet
    REF(1:length(n))=REF(1:length(n)).*0.5.*(1+cos(n)');                    %Fönstra den inspelade början
    
    clear REC m n  
else
    %Simulera dist
    disp('Simulating dist.')
    [b,a]=butter(3,[0.01 0.1]);
    DATA=y(:)+0.01*filter(b,a,y(:)).^2; %Simulera följande system
    REF=y(:);
    clear b a
end

%Zeropadda inspelningen
disp('Zeropadding')
DATA=[DATA' zeros(1,length(DATA))];
REF=[REF' zeros(1,length(REF))];
%Beräkna faltning av DATAkanalen
disp('Calculating convolution of DATA channel')
%Beräkna faltning av referenskanalen
out=fftfilt(yinv,DATA);
disp('Calculating convolution of reference channel')
ref=fftfilt(yinv,REF);
[apa, index0]=max(abs(ref)); %Hitta t=0 i referenskanalen

%Visa headroom
STRING=sprintf('Maximum value of DATA channel %f\nMaximum value of referens channel %f',max(abs(DATA)),max(abs(REF)));
disp(STRING)

%Beräkna fönster för varje delton

clear ref REF DATA apa
%grundtonen
L{1}=length(out)-index0;
if(L{1}>fs)
    L{1}=fs;
end

%Övertonerna
for n=2:MAXharmonic,
    L{n}=ceil(log(n)/time_k*fs);
end

INDEX{1}=index0-PREDELAY:index0-PREDELAY+L{1}-1;
INDEX{2}=index0-PREDELAY-L{2}:index0-PREDELAY-L{2}-1+round(0.5*L{2});
for n=3:MAXharmonic
    INDEX{n}=index0-PREDELAY-L{n}:index0-PREDELAY-L{n}-1+round(0.5*(L{n}-L{n-1}));
end

%Visa alla estimerade impulsvar

if(DISP_WINDOW==1)
    figure
end
for n=1:MAXharmonic
    part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
    Fpart{n}=fft(part{n});
    if(DISP_WINDOW==1)
        subplot(2,ceil(MAXharmonic/2),n),plot(0:1000/fs:1000*(length(part{n})-1)/fs,1000*part{n});
        STRING=sprintf('Imp. response; %d harmonic | %s',n,ELEMENT);
        title(STRING)
        xlabel('Time [ms]')
        grid on
    end    
end
if(DISP_WINDOW==1)
    drawnow
end


%Beräkna längsta möjliga fönster
MAXLENGTH=length(part{MAXharmonic});
if(ZOOM>MAXLENGTH)
    ZOOM=MAXLENGTH;
end


clear INDEX
%Visa inzoomad bild av grundtonens impulssvar
H_1=figure;plot(0:1000/fs:1000*(ZOOM-1)/fs,1000*part{1}(1:ZOOM));
title('IR-fundamental, Press mousebutton 4 times, (start-stop start-stop)');
xlabel('Time [ms]')
YLIM=ylim;
if(USE_WINDOW==0)
    N_WIND=zeros(4,1);
else
    N_WIND=ceil(fs/1000*WINDOW);
end

if(N_WIND(4)>MAXLENGTH)
    N_WIND(4)=MAXLENGTH;
    N_WIND(3)=MAXLENGTH-round(0.001*fs);
    disp('WARNING YOUR PREDEFINED WINDOW (WINDOW=[...]) can NOT be used')
    disp('WINDOW is too long')
    STRING=sprintf('The post window will start at %4.1f ms and end at %4.1f ms',N_WIND(4)/fs*1000,N_WIND(3)/fs*1000');
    disp(STRING)
    disp('Press ENTER to continue or Ctrl-C to halt')
    pause
end

n=1;

%Välj fönstringspunkter m h a musen
while( (0<=N_WIND(1) && N_WIND(1)<N_WIND(2)&& N_WIND(2)<N_WIND(3) && N_WIND(3)<N_WIND(4) && N_WIND(4)<=MAXLENGTH ) == 0 )
    disp('Press mousebutton 4 times')
    disp('First, press mousebutton at the beginning of the pre-window')
    disp('Second, press mousebutton at the end of the pre-window')
    disp('Third, press mousebutton at the beginning of the post-window')
    disp('Finally, press mousebutton at the end of the post-window')
    disp('T1<T2<T3<T4')
    disp('To achive a longer window, either the time T must be increased or the number MAXHarmonic must be decreased')
    
    if(n>1)
        text(1,0.90*YLIM(2),'ERROR - Unlegal window locations ! TRY AGAIN')
        disp('ERROR - Unlegal window locations ! TRY AGAIN')        
    end
    BUTTON=0;
    [Xg,Yg,BUTTON]=ginput(4);    %TRYCK 4 GGR MED MUSEN
    N_WIND=ceil(eps+Xg*fs/1000);
    n=n+1;
    if(min(N_WIND)<1)
        disp('Time must >= 0 TRY AGAIN')
    end
    if(max(N_WIND)>MAXLENGTH)
        disp('Time must > maxtime, TRY AGAIN')
    end
    
end



%Fönster till de valda fönstringspunkterna
WINDOW1=hann(2*(N_WIND(2)-N_WIND(1)))';
WINDOW2=hann(2*(N_WIND(4)-N_WIND(3)))';

%Applicera ovan fönster på alla deltoner
for n=1:MAXharmonic
    TRUNK{n}=part{n}(N_WIND(1):N_WIND(4));
    TRUNK{n}(1:N_WIND(2)-N_WIND(1))=TRUNK{n}(1:N_WIND(2)-N_WIND(1)).*WINDOW1(1:N_WIND(2)-N_WIND(1));
    TRUNK{n}(N_WIND(3)-N_WIND(1):N_WIND(4)-N_WIND(1))=TRUNK{n}(N_WIND(3)-N_WIND(1):N_WIND(4)-N_WIND(1)).*WINDOW2(N_WIND(4)-N_WIND(3):end);
    LENGTH=2*ceil(length(TRUNK{n})/2);
    Ftrunk{n}=fft(TRUNK{n},LENGTH); %Skapa FFT med jämn längd, viktigt!
end

%Plotta den trunkerade grundtonens impulsvar (innan första reflexen) +fönstrena
if(DISP_WINDOW==1)
    MAX=max(abs(TRUNK{1}));
    figure,plot(0:1000/fs:1000*(length(TRUNK{1})-1)/fs,TRUNK{1},...
        0:1000/fs:1000*(length(WINDOW1)/2-1)/fs,MAX.*WINDOW1(1:length(WINDOW1)/2),'r:',...
        N_WIND(3)*1000/fs-N_WIND(1)*1000/fs:1000/fs:N_WIND(3)*1000/fs-N_WIND(1)*1000/fs+1000*(length(WINDOW2)/2)/fs,MAX.*WINDOW2(length(WINDOW2)/2:end),'r:' )
    legend('Impulse','Raised cosine windows')
    xlabel('Time [ms]')
    title('Truncated impulse response')
    %Spara plot till fil
    STRING=sprintf('%s%s_IMP.bmp',PATH,ELEMENT);
    disp(['Saving plot to file: ' STRING]);
    print('-dbmp16m',STRING) 
    
    %Plotta en Energy Time Curve för hela resultatet, samt markera var
    %deltonerna kan hittas
    t=-index0/fs:1/fs:(length(out)-index0-1)/fs ;
    figure, plot(t,10*log10(out.^2+eps));
    title('Energy Time Curve');
    YLIM=ylim;
    ylim([YLIM(2)-100 YLIM(2)]);
    grid on
    xlabel('Time [s]');
    ylabel('Level [dB]');
    hold on
    plot([0 0],[YLIM(2)-100 YLIM(2)],'k','linewidth',2);
    
    for i=2:MAXharmonic
        plot([-L{i}/fs -L{i}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
    end
    
    hold off
    legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic');
    drawnow
end

clear L t out

%Beräkna frekvenserna för alla deltonerna, de trunkera deltonerna, samt
%beräkna nivån i frekvensdomänen för de trunkerade impulsvaren för varje
%delton
for n=1:MAXharmonic
    F{n}=linspace(0,fs/(2*n),length(Fpart{n})/2);
    FT{n}=linspace(0,fs/2/n,length(Ftrunk{n})/2);
    dBT{n}=20*log10(abs(Ftrunk{n}(1:length(Ftrunk{n})/2)));
end

%Skapa konstant upplösning i frekvensled på logskala
ktemp=log(f_end/f_start);
Ftemp=f_start*exp(ktemp*linspace(0,1,2*RES));
COUNT=ones(MAXharmonic,1);
POWER=zeros(2*RES,1);
k1=1;
%Sortera in energin från linjärskala till logskala
for m=2:2:length(Ftemp)-1
    
    for n=1:MAXharmonic
        I=find(and(F{n}>Ftemp(m-1), F{n}<Ftemp(m+1))); 
        if(isempty(I)~=1)
            dB{n}(COUNT(n))=10*log10(mean(abs(Fpart{n}(I)).^2));
            Flog{n}(COUNT(n))=Ftemp(m);
            COUNT(n)=COUNT(n)+1;
            if(n>1)
                POWER(k1)=POWER(k1)+mean(abs(Fpart{n}(I)).^2);
                Fpower(k1)=Ftemp(m);
            end
            
        end
    end
    k1=k1+1;
end

%Kalibreringssekvens
if(Kalibrera==1)
    [apa, I]=find(Flog{1}>CALFRERQ);
    disp('KALIBRERING, MATA IN NEDANSTÅENDE VÄRDE TILL PARAMETERN OFFSET_dB')
    disp(dB{1}(I(1))) % Kalibrering
    return  %Skriv in värdet på OFFSET_dB
else
    CAL=CALdB-OFFSET_dB;
end


clear ktemp Ftemp COUNT I

%Plotta frekvensgången för grundton + 2:a & 3:e delton
figure,semilogx(...
    Flog{1},dB{1}+CAL,'k:',...
    FT{1},dBT{1}+CAL,'k',...
    Flog{2},dB{2}+CAL,'r:',...
    FT{2},dBT{2}+CAL,'r',...
    Flog{3},dB{3}+CAL,'b:',...
    FT{3},dBT{3}+CAL,'b','linewidth',2);
%Fpower,10*log10(POWER(1:length(Fpower))+eps)+CAL,'g:',

grid on


ylim([20 100]);
if(TICK==1)
    set(gca,'XTick',XTICK);
    set(gca,'XTickLabel',XTICK);
    xlim([XTICK(1) XTICK(end)]);
else
    xlim([f_start f_end]);
end



STRING=sprintf('Frequency [Hz] | Trolldist v. %2.2f Beta',Trollver);
%STRINGY=sprintf('SPL [dB] @ %d cm ref 20 uPa, %s',MIKAVSTAND,VOLTAGE);
STRINGY=sprintf('SPL [dB] @ %d cm !UNCALIBRATED!, %s',MIKAVSTAND,VOLTAGE);
xlabel(STRING);
ylabel(STRINGY);
legend('Fundamental','Truncated before first reflection','2:nd harmonic',...
    'Truncated before first reflection','3:rd harmonic','Truncated before first reflection',0);
STRING=sprintf('Dist. and Freq resp. %s | %s',CONDITION,ELEMENT);
title(STRING)
set(gcf,'Position',[100 100 480 480]);
drawnow

%Spara plot till fil
STRING=sprintf('%s%s.bmp',PATH,ELEMENT);
disp(['Saving plot to file: ' STRING]);
print('-dbmp16m',STRING) 


%SPARA impulsvaren till disk
IR=TRUNK;
time=0:1/fs:(length(IR{1})-1)/fs;
disp(['Saving impulse responses to file' STRING]);
save(ELEMENT,'IR','fs','time')

%Beräkna CDS Cumulative Decay Spectrum
[Bs,Fs,Ts]=specgram(IR{1}*10^(CAL/20),64,fs,hann(64),60);
figure,waterfall(-Fs,1000*Ts,20*log10(abs(Bs')+10^(30/20)))
xlim([-20000 0])
xlabel('Frekvens [Hz]')
ylabel('Tid [s]')
zlabel('SPL @ 50 cm [dB]')
view(160,30)
zlim([30 100])
grid on
title('Cumulative Spectral Deyac Prototype')
XTICK=[0;5000;10000;15000;20000];
set(gca,'XTick',-XTICK(end:-1:1));
set(gca,'XTickLabel',XTICK(end:-1:1));
drawnow

if(TESTcrossover==1)
    
    Length1=1024;   %Längden på minimization of sum |B-H*A|^2*Wt nr1
    Length2=256;    %%Längden på minimization of sum |B-H*A|^2*Wt nr1
    Length_Hilbert=4096;    %Längden på Hilbert transform kompensatorn OBS bör gärna vara minst 2*length(IR{1}) OBS måste vara större eller lika stor som Length_levinson!!!
    Length_ARyulewalk=1024; 
    Length_levinson=256;    %Längden på Levinson-Durbin kompensatorn (med modifierad akf)
    
    FFT=fft(IR{1},Length_Hilbert);  %Används till Hilbert transformen
    
    [Hcomplex1,Wcomplex1]=freqz(IR{1},1,Length1);
    [Hcomplex2,Wcomplex2]=freqz(IR{1},1,Length2);
    FH=linspace(0,fs/2,Length_Hilbert/2);
    [APA,Ih]=find(FH>Comp_high); %Övre gränsfrekvens för Hilbert kompensatorn
    [APA,Il]=find(FH<Comp_low); %Undre gränsfrekvens för Hilbe kompensatorn
    
    %Skapa den inverterade frekvensgången som ska användas till
    %Hilberttransformen
    Finv=1./FFT(1:Length_Hilbert/2);
    Finv(Il)=Finv(Il(end)); 
    Finv(Ih)=Finv(Ih(1));
    
    FAR=FFT;
    FAR(Il)=FAR(Il(end)); 
    FAR(Ih)=FAR(Ih(1));
    R=FAR.*conj(FAR);   %Akf i frekvensdomänen
    r=real(ifft(R));    %Akf i tidsdomänen
    disp('Calculating AR model from autocorrelation sequence based on Levinson-Durbin Recursion')
    EQminphase_levinson=levinson(r,Length_levinson-1);
    Hlevinson=freqz(conv(IR{1},EQminphase_levinson),1,64);
    EQminphase_levinson=EQminphase_levinson./median(abs(Hlevinson));
    
    
    %     [H,W]=freqz(IR{1},1,1024); 
    %     [Index,APA]=find(W/pi>18000/(0.5*fs)); 
    %     H(Index)=H(Index(1)); 
    %     [Index,APA]=find(W/pi<200/(0.5*fs)); 
    %     H(Index)=H(Index(end)); 
    %     [B,A]=invfreqz(H,W,0,1024); 
    %     h_inv=A/B;
    %     
    
    [Ih,APA]=find(Wcomplex2/pi>Comp_high/(0.5*fs));  %Övre gränsfrekvens för kompensatorn
    [Il,APA]=find(Wcomplex2/pi<Comp_low/(0.5*fs)); %Undre gränsfrekvens för komensatorn
    Hcomplex2(Il)=Hcomplex2(Il(end)+1); 
    Hcomplex2(Ih)=Hcomplex2(Ih(1)-1);
    
    [Ih,APA]=find(Wcomplex1/pi>Comp_high/(0.5*fs));  %Övre gränsfrekvens för kompensatorn
    [Il,APA]=find(Wcomplex1/pi<Comp_low/(0.5*fs)); %Undre gränsfrekvens för komensatorn
    Hcomplex1(Il)=Hcomplex1(Il(end)); 
    Hcomplex1(Ih)=Hcomplex1(Ih(1));
    
    %h_new=real(ifft([Finv Finv(end:-1:1)]));
    %h_new=[h_new(:)' zeros(1,length(h_new))];
    %Interpolera värden om man vill
    %Finv=interp(Finv,8);
    
    % :) Skapa impulsvaret för det minimum-fas system som innehar den
    % inverterade frekvensgången baserat på Hilbert transform  :) 
    %Lilltrolls favoritrad i MATLAB :)
    disp('Calculating inverse based on Hilbert Transform')
    EQminphase_hilbert=real(ifft(exp(conj(hilbert(log(abs([Finv Finv(end:-1:1)])))))));
    Hhilbert=freqz(conv(IR{1},EQminphase_hilbert),1,64);
    EQminphase_hilbert=EQminphase_hilbert./median(abs(Hhilbert)); %Justera så att ampltuden blir 1 hos systemet faltat med inversen
    EQminphase_hilbert=EQminphase_hilbert(1:Length_Hilbert);
    %Skapa impulsvaret för det minimum-fas system som innehar den
    %inverterade frekvensgången baserat på AR Yule-Walker
    
    disp('Calculating inverse based on AR Yule-Walker')
    disp('This can take some time')
    EQminphase_ARyulewalk=aryule(IR{1},length(IR{1}));
    Haryulewalk=freqz(conv(IR{1},EQminphase_ARyulewalk),1,64);
    EQminphase_ARyulewalk=EQminphase_ARyulewalk./median(abs(Haryulewalk)); %Justera så att ampltuden blir 1 hos systemet faltat med inversen
    
    %Skapa impulsvaret för det minimum-fas system som innehar den
    %inverterade frekvensgången baserat på invfreqz
    
    disp('Calculating inverse based on minimization of sum |B-H*A|^2*Wt over the freqs in W.')
    disp('This might take even longer time and a lot of memory !!!')
    
    
    Vikt1=ones(1,Length1);
    [APA,EQminphase_invfreqz1]=invfreqz(Hcomplex1,Wcomplex1,0,Length1-1,Vikt1);
    Hinvfreqz1=freqz(conv(IR{1},EQminphase_invfreqz1),1,64);
    EQminphase_invfreqz1=EQminphase_invfreqz1./median(abs(Hinvfreqz1)); %Justera så att ampltuden blir 1 hos systemet faltat med inversen
    [Btest,Atest]=invfreqz(Hcomplex1,Wcomplex1,4,4,Vikt1);
    
    Vikt2=ones(1,Length2);
    Vikt2(1)=0.4;   %Handpåläggning för att vikta basen, typiskt <1
    Vikt2(2)=0.5;   %Handpåläggning för att vikta basen, typiskt <1
    [APA,EQminphase_invfreqz2]=invfreqz(Hcomplex2,Wcomplex2,0,Length2-1,Vikt2);
    Hinvfreqz2=freqz(conv(IR{1},EQminphase_invfreqz2),1,64);
    EQminphase_invfreqz2=EQminphase_invfreqz2./median(abs(Hinvfreqz2)); %Justera så att ampltuden blir 1 hos systemet faltat med inversen
    
    figure,plot(0:1000/fs:1000*(length(EQminphase_hilbert)-1)/fs,EQminphase_hilbert,'o',0:1000/fs:1000*(length(EQminphase_ARyulewalk)-1)/fs,EQminphase_ARyulewalk,'o',0:1000/fs:1000*(length(EQminphase_invfreqz1)-1)/fs,EQminphase_invfreqz1,'o',0:1000/fs:1000*(length(EQminphase_invfreqz2)-1)/fs,EQminphase_invfreqz2,'o',0:1000/fs:1000*(length(EQminphase_levinson)-1)/fs,EQminphase_levinson,'o')
    xlabel('Time [ms]');
    title('Minimum-phase inverses');
    STRING1=sprintf('Based on minimization of sum |B-H*A|^2*Wt with length %d',Length1);
    STRING2=sprintf('Based on minimization of sum |B-H*A|^2*Wt with length %d',Length2);
    STRING3=sprintf('Based on Levinson-Durbin Recursion from modified acf with length %d',Length_levinson);
    STRING4=sprintf('Based on Hilbert transform %d',Length_Hilbert);
    legend(STRING4,'Based on AR Yule-Walker',STRING1,STRING2,STRING3);
    disp('Saving FIR-equlizer coeff. to disk');
    save(['EQminphase-' ELEMENT],'fs','EQminphase_hilbert','EQminphase_ARyulewalk','EQminphase_invfreqz1','EQminphase_invfreqz2') 
    clear eit h_aryule
    drawnow
    %---------------------------
    %Labb område för att plotta prestandan av inverserna
    %---------------------------
    Moreplots=1;
    if(Moreplots==1)
        RESULT_hilbert=conv(TRUNK{1},EQminphase_hilbert); %Resultatet av impulsvaret faltat med digital-EQ baserad på Hilbert
        RESULT_levinson=conv(TRUNK{1},EQminphase_levinson); %%Resultatet av impulsvaret faltat med digital-EQ baserad på AR Yule-Walker
        RESULT_invfreqz1=conv(TRUNK{1},EQminphase_invfreqz1); %%Resultatet av impulsvaret faltat med digital-EQ baserad på invfreqz
        RESULT_invfreqz2=conv(TRUNK{1},EQminphase_invfreqz2);
        
        figure,plot(0:1000/fs:1000*(length(RESULT_hilbert)-1)/fs,RESULT_hilbert,0:1000/fs:1000*(length(RESULT_levinson)-1)/fs,RESULT_levinson,0:1000/fs:1000*(length(RESULT_invfreqz1)-1)/fs,RESULT_invfreqz1,0:1000/fs:1000*(length(RESULT_invfreqz2)-1)/fs,RESULT_invfreqz2);
        xlabel('Time [ms]');
        title('Minimum-phase invers + Impulse response');
        legend(STRING4,STRING3,STRING1,STRING2);
        drawnow
        [H_RESULT_Hilbert,W_RESULT]=freqz(RESULT_hilbert,1,Length1,fs);
        [H_RESULT_levinson,W_RESULT]=freqz(RESULT_levinson,1,Length1,fs);
        [H_RESULT_invfreqz1,W_RESULT]=freqz(RESULT_invfreqz1,1,Length1,fs);
        [H_RESULT_invfreqz2,W_RESULT]=freqz(RESULT_invfreqz2,1,Length1,fs);        
        
        figure,semilogx(W_RESULT,20*log10(abs(H_RESULT_Hilbert)),W_RESULT,20*log10(abs(H_RESULT_levinson)),W_RESULT,20*log10(abs(H_RESULT_invfreqz1)),W_RESULT,20*log10(abs(H_RESULT_invfreqz2)),'linewidth',2)
        line(Comp_low*[1 1],ylim,'linewidth',3,'color',[0 0 0]);
        line(Comp_high*[1 1],ylim,'linewidth',3,'color',[0 0 0]);
        Title('Frequency response performance of compensator')
        xlabel('Frequency [Hz]');
        ylabel('Level [dB]');
        legend(STRING4,STRING3,STRING1,STRING2);
        grid on
        axis tight
        ylim([-5 5]);
        drawnow
        disp('Saving FIR koeff for DSP')
        SAVEDATA=EQminphase_levinson(end:-1:1)/max(abs(EQminphase_levinson));
        save EQlevinson.txt SAVEDATA -tabs -ascii
        
        % EJ KLART   
        %     TOTAL=conv(eit,RESULT(end:-1:1));
        %     [APA,INDEXc]=max(abs(TOTAL));
        %     LENGTH=512;
        %     WINDOW=hann(LENGTH);
        %     EQ=TOTAL(INDEXc-LENGTH/2+1:INDEXc+LENGTH/2).*WINDOW';
        %         
        %     figure,plot(0:1000/fs:1000*(length(EQ)-1)/fs,EQ)
        %     xlabel('Time [ms]');
        %     title('Minimum-phase + Excess-phase invers');
        %     save(['EQ-' ELEMENT],'fs','EQminphase','EQ') 
        %     
        %     figure,plot(0:1000/fs:1000*(length(RESULT)-1)/fs,RESULT)
        %     xlabel('Time [ms]');
        %     title('Minimum-phase invers + Impulse response');
        %     %save(['EQ' ELEMENT],'fs','EQminphase') 
        % 
        %     [B,A]=butter(3,2200/(0.5*fs))
        %     TOTALcrossover=filter(B,A,EQ);
        %     figure,plot(0:1000/fs:1000*(length(TOTALcrossover)-1)/fs,TOTALcrossover)
        %     xlabel('Time [ms]');
        %     title('Minimum-phase + Excess-phase invers + Crossover [3:th order Butterworth @ 2.2kHz]');
        %     
        %     TOTALcrossoverMP=filter(B,A,EQminphase);
        %     figure,plot(0:1000/fs:1000*(length(TOTALcrossoverMP)-1)/fs,TOTALcrossoverMP)
        %     xlabel('Time [ms]');
        %     title('Minimum-phase EQ + Crossover [3:th order Butterworth @ 2.2kHz]');
        %     
        %     %Beräkna resultatet
        %     FINALRESULT_MP=conv(TOTALcrossoverMP,TRUNK{1});
        %     FINALRESULT=conv(TOTALcrossover,TRUNK{1});
        %     Time=0:1000/fs:1000*(length(FINALRESULT)-1)/fs;
        %     Time_MP=0:1000/fs:1000*(length(FINALRESULT_MP)-1)/fs;
        %     figure,plot(Time_MP,FINALRESULT_MP)
        %     title('Result of digital crossover with minimumphase EQ')
        %     xlabel('Time [ms]');
        %     figure,freqz(FINALRESULT_MP,1,length(FINALRESULT),fs)
        %     title('Result of digital crossover with minimumphase EQ')
        %     set(gca,'xscale','log')
        %     ylim([-40 10])
        %     xlim([100 20000])
        %     
        %     figure,plot(Time,FINALRESULT)
        %     title('Result of digital crossover with EQ')
        %     xlabel('Time [ms]');
        %     figure,freqz(FINALRESULT,1,length(FINALRESULT),fs)
        %     title('Result of digital crossover with EQ')
        %     set(gca,'xscale','log')
        %     ylim([-40 10])
        %     xlim([100 20000])
        %     
        %     SAVEDATA=EQminphase/max(abs(EQminphase));
        %     
        %     save EQ4.txt SAVEDATA -tabs -ascii
        %     save EQB.txt B -tabs -ascii
        %     save EQA.txt A -tabs -ascii
        %     
        %fid = fopen('EQ2.txt','wb')
        %fprintf(fid,'%.10f\t',RESULT/norm(RESULT,inf));
        %fclose(fid)
        
        %----------------------------------        
        %Testområde för olika inversmetoder
        %----------------------------------
    end
    
    TEST=0;
    if(TEST==1)
        h_hilbert=real(ifft(exp(conj(hilbert(log(abs(1./fft(IR{1}))))))));
        L=length(IR{1});
        h_aryule=aryule(IR{1},L);
        K=mean(abs(freqz(conv(IR{1},h_aryule),1,256,fs)));
        h_arburg=arburg(IR{1},L-1);
        h_arcov=arcov(IR{1},floor(L/2));
        [H,W]=freqz(IR{1},1,1024);
        [Index,APA]=find(W/pi>18000/(0.5*fs));
        H(Index)=H(Index(1));
        [Index,APA]=find(W/pi<200/(0.5*fs));
        H(Index)=H(Index(end));
        [B,A]=invfreqz(H,W,0,1024);
        h_inv=A/B;
        
        plot(0:1000/fs:1000*(length(h_hilbert)-1)/fs,h_hilbert...
            ,0:1000/fs:1000*(length(h_aryule)-1)/fs,h_aryule/K...
            ,0:1000/fs:1000*(length(h_arburg)-1)/fs,h_arburg...
            ,0:1000/fs:1000*(length(h_arcov)-1)/fs,h_arcov);
        legend('Hilbert Transform','AR parameter estimation via Yule-Walker method'...
            ,'AR parameter estimation via Burg method.','AR parameter estimation via covariance method')
        xlabel('Tid [ms]')
        
        [H_hilbert,W_hilbert]=freqz(conv(IR{1},h_hilbert),1,512,fs);
        [H_aryule,W_aryule]=freqz(conv(IR{1},h_aryule),1,512,fs);
        [H_inv,W_inv]=freqz(conv(IR{1},h_inv),1,512,fs);
        
        figure,semilogx(W_hilbert,20*log10(abs(H_hilbert)),W_aryule,20*log10(abs(H_aryule/K)),W_inv,20*log10(abs(H_inv)));
        grid on
        
        semilogx(W_inv,20*log10(abs(H_inv)));
        grid on
    end
    
    %Använd Prony för att skapa en ARMA modell utifrån minfas MA modellen.
    USEPRONY=0;
    if(USEPRONY==1)
        clear H W SIZE
        for k=1:4
            SIZE(k)=2^(k+1);
            [B,A]=prony(eit,SIZE(k),SIZE(k));
            [H(:,k),W(:,k)]=freqz(conv(IR{1},B),A,fs,fs);    
        end
        figure,semilogx(W,20*log10(abs(H)),'linewidth',2);
        grid on
        xlabel('Frekvens [Hz]');
        ylabel('Nivå [dB]');
        legend(num2str(SIZE'))
        xlim([200 10000])
        ylim([-5 5])
        title('Prony´s metod') 
        figure, zplane(B,A)
        title('Prony´s metod') 
        
        %----------------------------------------
        %Beräkna Zrummet som 3D plot
        %----------------------------------------
        [R,I]=meshgrid(-1:0.003:1,-1:0.003:1); %Skapa ett komplex talplan
        C=R+j*I;
        [xx,yy]=find(abs(C)>1); %Hitta alla punkter innanför enhetcirkeln
        w=linspace(0,2*pi,360);
        X=0;
        R=0;
        RAND=exp(j*w); %enhetscirkeln
        Brand = polyval(B,RAND);    %Beräkna tälarpolynomet på enhetcirkeln
        Arand= polyval(A,RAND);     %Beräkna nämnarpolynomet på enhetcirkeln
        MagRAND=20*log10((abs(Brand./Arand))); %Beräkna frekvensgången i dB (på enhetscirkeln)
        Broot=roots(B); %Hitta rötterna till täljarpolynomet
        Aroot=roots(A); %Hitta rötterna till nämnarpolynomet
        Bc = polyval(B,C);  %Beräkna tälarpolynomet innanför enhetcirkeln
        Ac= polyval(A,C);   %Beräkna nämnarpolynomet innanför enhetcirkeln
        Mag=20*log10((abs(Bc./Ac)));    %Beräkna tältduken uttryckt i dB
        for i=1:length(xx)
            Mag(xx(i),yy(i))=NaN;   %Definiera alla punkter utanför enhetcirkeln som NaN
        end
        %plotta duken
        figure,surf(R,I,Mag,'FaceColor','interp',... 
            'EdgeColor','none',...
            'FaceLighting','phong')
        view(135,30);   %Sätt kameravinkel
        camlight(-30,15)    %Sätt ljusvinkeln
        hold on
        plot3(real(RAND),imag(RAND),MagRAND,'r','linewidth',3); %Plotta nivåkurvan på enhetscirkeln
        plot3(real(RAND),imag(RAND),zeros(length(RAND)),'k','linewidth',3) %Plotta enhetscirkeln
        plot(real(Broot),imag(Broot),'or')  %Plotta nollställena
        plot(real(Aroot),imag(Aroot),'xr')  %Plotta polerna
        hold off
        lighting gouraud
        material shiny
        shading interp
        
        xlabel('Real');
        ylabel('Imag');
        title('Z-plane');
        zlabel('Level [dB]');
        zlim([0 100])
        grid on
        text(-1,0,5,sprintf('%d Hz',fs/2),'BackgroundColor',[1 1 1],'Color',[1 0 1]);
        text(0,-1,5,sprintf('-%d Hz',fs/4),'BackgroundColor',[1 1 1],'Color',[1 0 1]);
        text(0,1,5,sprintf('%d Hz',fs/4),'BackgroundColor',[1 1 1],'Color',[1 0 1]);
        text(1,0,5,sprintf('0 Hz',fs/4),'BackgroundColor',[1 1 1],'Color',[1 0 1]);
        
    end
    %Testar digitala delningsfilter + digital EQ. Kräver DAQ toolboxen +
    %externa m filer
    REALTIME_EQ=0;  %OBS fungerar inte utan funktionerna qmoredataFFT qmoredataIIR, fråga Lilltroll efter dessa
    if(REALTIME_EQ==1)
        
        METHOD='FIR_FFT' %'IIR' eller 'EJ KLAR FIR_FFT+IIR EJ KLAR' %Välj filtreringstyp
        
        global Zi A B % Tillfälligt fulhack, enkelt med lite fult
        fs=44100;
        ai = analoginput('winsound');
        addchannel(ai,1);
        ao = analogoutput('winsound');
        addchannel(ao,1);
        %You should use manual triggers when starting multiple device objects because this trigger type executes faster than other trigger types with the exception of hardware triggers. Additionally, to synchronize the input and output of data, you should configure the ManualTriggerHwOn property to Trigger for ai. 
        set([ai ao],'TriggerType','Manual')
        set(ai,'ManualTriggerHwOn','Trigger')
        %Configure ai for continuous acquisition, call the callback function qmoredata whenever 1000 samples are output, and call daqcallback when ai and ao stop running.
        if(METHOD=='FIR_FFT')
            FFTsize=4096*2; %Längden på varje datapaket, bör vara en multipel på 2 samt minst 2ggr större än Length_FIR
            Length_FIR=1024; %Längden på 
            B=fft(0.15*EQminphase_invfreqz1(1:Length_FIR),FFTsize);
            A=1;
            Zi=zeros(1,Length_FIR);
            set(ao,'SamplesOutputFcn',{'qmoredataFFT',ai,FFTsize,Length_FIR});    
        elseif(METHOD=='IIR')
            [B,A]=butter(3,0.25);   %Testfilter
            set(ao,'SamplesOutputFcn',{'qmoredataIIR',ai});    
        end
        
        set([ai ao],'SampleRate',fs);
        fs_real=get(ai,'SampleRate');
        disp(sprintf('Verklig fs=%d Hz',fs_real))
        set(ai,'SamplesPerTrigger',Inf);
        set(ao,'SamplesOutputFcnCount',FFTsize-Length_FIR);
        set([ai ao],'StopFcn',@daqcallback);
        
        %Queue data in the engine, start the device objects, and execute the manual triggers. data = zeros(4000,1);
        data = zeros(fs,1);
        putdata(ao,data);
        start([ai ao]);
        trigger([ai ao]);
        
        stop([ai ao])
        %-------------------------------
        
    end
    
end
    
    %DISTMÄTNING VID KONSTANT UTNIVÅ
    if(Konstant_level>0)
        
        %Hitta frekvensområdet som ska nivåkompenseras
        [APA,Il]=find(FT{1}<Konstant_level_start); %Hz
        [APA,Ih]=find(FT{1}>Konstant_level_end); %Hz
        
        %Skapa magnituden på frekvensgången från de trunkerade impulsvaren
        Fmag=abs(Ftrunk{1}./Ftrunk{1}(Il(end)));
        Fmag=Fmag(1:length(Fmag)/2);
        
        %Skapa den inverterade frekvensgången
        Finv=1./Fmag;
        Finv(Il)=linspace(0.1,1,length(Il)); %Minska amplituden linjärt innan området börjar
        Finv(Ih)=Finv(Ih(1));
        
        %Interpolera värden om man vill
        %Finv=interp(Finv,8);
        
        % :) Skapa impulsvaret för det minimum-fas system som innehar den
        % inverterade frekvensgången baserat på Hilbert transform  :) 
        %Lilltrolls favoritrad i MATLAB :)
        eit=real(ifft(exp(conj(hilbert(log([Finv Finv(end:-1:1)]))))));
        %Spara värderna som kan användas för digital minimum-fas EQ
        STRING=sprintf('%s_EQ',ELEMENT);
        disp(['Saving EQ coef. to file' STRING]);
        save(STRING,'eit');
        
        %Beräkna det minimum-fas EQ ade chirpen
        yB=fftfilt(eit,y);
        clear FPOWER
        
        %Korrigera till önskad nivå
        for LOOP=1:length(Konstant_level_LEVEL)
            
            AMPdB=Konstant_level_LEVEL(LOOP)-20*log10(abs(Ftrunk{1}(Il(end))))-CAL;
            STRING=sprintf('Testing dist. at %3.1f dB',Konstant_level_LEVEL(LOOP));
            disp(STRING)
            PLAY=10.^(AMPdB/20)*0.99*AMP*yB;
            
            %Kontrollera att signalen inte klipper
            MaxA=max(abs(PLAY));
            if(MaxA>0.99)
                disp('Output volume is not high enougth to achivethe wanted SPL')
                disp(MaxA)
                LOOP=LOOP-1;
                return
            end
            
            %Spela in signalen
            disp('Recording signal')
            pause(0.1);
            wavplay([PLAY(:) AMP*y(:)],fs,'async');REC=wavrecord(length(PLAY),fs,2);
            disp('Processing recording');
            
            
            if(max(max(abs(REC)>0.99)))
                disp('!WARNING!, input signal is overloaded !!!Program halted!!!')
                break
            end
            
            if(max(max(abs(REC)>0.8)))
                disp('!WARNING!, input signal is probably overloaded')
                disp('Consider lower the gain to the soundcard input/Mixer levels (SYSTEM MUST THEN BE RECALIBRATED)')
                disp('Press ENTER to continue anyway, Press Crtl-C to stop')
                pause
            end
            
            
            DATA=REC(:,1).*tukeywin(length(REC),0.2);
            REF=REC(:,2).*tukeywin(length(REC),0.2);
            DATA=DATA(:).*tukeywin(length(DATA),0.2);
            REF=REF(:).*tukeywin(length(REF),0.2);
            clear REC
            
            disp('Zeropadding')
            DATA=[DATA' zeros(1,length(DATA))];
            REF=[REF' zeros(1,length(REF))];
            disp('Calculating convolution of DATA channel')
            out=fftfilt(yinv,DATA);
            disp('Calculating convolution of reference channel')
            ref=fftfilt(yinv,REF);
            
            [apa, i]=max(abs(ref)); %Jämför med referenssignalen
            
            STRING=sprintf('Maximum value of DATA channel %f\nMaximum value of referens channel %f',max(abs(DATA)),max(abs(REF)));
            disp(STRING)
            
            %grundtonen
            L{1}=length(out)-i;
            if(L{1}>fs)
                L{1}=fs;
            end
            for n=2:MAXharmonic,
                L{n}=ceil(log(n)/time_k*fs);
            end
            
            INDEX{1}=i-PREDELAY:i-PREDELAY+L{1}-1;
            INDEX{2}=i-PREDELAY-L{2}:i-PREDELAY-L{2}-1+round(0.5*L{2});
            for n=3:MAXharmonic
                INDEX{n}=i-PREDELAY-L{n}:i-PREDELAY-L{n}-1+round(0.5*(L{n}-L{n-1}));
            end
            
            if(DISP_WINDOW==1)
                figure
            end
            for n=1:MAXharmonic
                part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
                Fpart{n}=fft(part{n});
                if(DISP_WINDOW==1)
                    subplot(2,ceil(MAXharmonic/2),n),plot(0:1000/fs:1000*(length(part{n})-1)/fs,1000*part{n});
                    STRING=sprintf('Imp. response; %d harmonic',n);
                    title(STRING)
                    xlabel('Time [ms]')
                end    
            end
            
            if(DISP_WINDOW==1)
                drawnow
            end
            
            
            clear INDEX
            
            
            H_1=figure;plot(0:1000/fs:1000*(ZOOM-1)/fs,1000*part{1}(1:ZOOM));
            title('IR-fundamental, Press mousebutton 4 times, (start-stop start-stop)');
            xlabel('Time [ms]')
            YLIM=ylim;
            drawnow
            n=1;
            while( ( N_WIND(1)<N_WIND(2)&& N_WIND(2)<N_WIND(3) && N_WIND(3)<N_WIND(4) ) == 0 )
                disp('Press mousebutton 4 times')
                disp('First, press mousebutton at the beginning of the pre-window')
                disp('Second, press mousebutton at the end of the pre-window')
                disp('Third, press mousebutton at the beginning of the post-window')
                disp('Finally, press mousebutton at the end of the post-window')
                disp('T1<T2<T3<T4')
                
                if(n>1)
                    text(1,0.90*YLIM(2),'ERROR - Unlegal window locations ! TRY AGAIN')
                    disp('ERROR - Unlegal window locations ! TRY AGAIN')        
                end
                BUTTON=0;
                [Xg,Yg,BUTTON]=ginput(4);    %TRYCK 4 GGR MED MUSEN
                N_WIND=round(Xg*fs/1000);
                n=n+1;
            end
            
            WINDOW1=hann(2*(N_WIND(2)-N_WIND(1)))';
            WINDOW2=hann(2*(N_WIND(4)-N_WIND(3)))';
            
            for n=1:MAXharmonic
                TRUNK{n}=part{n}(N_WIND(1):N_WIND(4));
                TRUNK{n}(1:N_WIND(2)-N_WIND(1))=TRUNK{n}(1:N_WIND(2)-N_WIND(1)).*WINDOW1(1:N_WIND(2)-N_WIND(1));
                TRUNK{n}(N_WIND(3)-N_WIND(1):N_WIND(4)-N_WIND(1))=TRUNK{n}(N_WIND(3)-N_WIND(1):N_WIND(4)-N_WIND(1)).*WINDOW2(N_WIND(4)-N_WIND(3):end);
                FtrunkDIST{n}=fft(TRUNK{n});
            end
            MAX=max(abs(TRUNK{1}));
            if(DISP_WINDOW==1)
                
                figure,plot(0:1000/fs:1000*(length(TRUNK{1})-1)/fs,TRUNK{1},...
                    0:1000/fs:1000*(length(WINDOW1)/2-1)/fs,MAX.*WINDOW1(1:length(WINDOW1)/2),'r:',...
                    N_WIND(3)*1000/fs-N_WIND(1)*1000/fs:1000/fs:N_WIND(3)*1000/fs-N_WIND(1)*1000/fs+1000*(length(WINDOW2)/2)/fs,MAX.*WINDOW2(length(WINDOW2)/2:end),'r:' )
                legend('Impulse','Raised cosine windows')
                xlabel('Time [ms]')
                title('Truncated impulse response')
                
                t=-index0/fs:1/fs:(length(out)-index0-1)/fs ;
                figure, plot(t,10*log10(out.^2+eps));
                title('Energy Time Curve');
                YLIM=ylim;
                ylim([YLIM(2)-100 YLIM(2)]);
                grid on
                xlabel('Time [s]');
                ylabel('Level [dB]');
                hold on
                plot([0 0],[YLIM(2)-100 YLIM(2)],'k','linewidth',2);
                
                for i=2:MAXharmonic
                    plot([-L{i}/fs -L{i}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
                end
                
                hold off
                legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic');
                drawnow
            end
            clear L t out
            
            for n=1:MAXharmonic
                F{n}=linspace(0,fs/(2*n),length(Fpart{n})/2);
                FT{n}=linspace(0,fs/2/n,length(FtrunkDIST{n})/2);
                dBT{n}=20*log10(abs(FtrunkDIST{n}(1:length(FtrunkDIST{n})/2)));
            end
            
            
            [apa, I]=find(F{1}>CALFRERQ);
            clear apa Fpower
            %CAL=CALdB-20*log10(abs(Fpart{1}(I(1)))); % Kalibrering
            
            ktemp=log(f_end/f_start);
            Ftemp=f_start*exp(ktemp*linspace(0,1,2*RES));
            COUNT=ones(MAXharmonic,1);
            POWER=zeros(length(FT{2}),1);
            k1=1;
            for m=2:2:length(Ftemp)-1
                
                for n=1:MAXharmonic
                    I=find(and(F{n}>Ftemp(m-1), F{n}<Ftemp(m+1))); 
                    if(isempty(I)~=1)
                        dB{n}(COUNT(n))=10*log10(mean(abs(Fpart{n}(I)+eps).^2));
                        Flog{n}(COUNT(n))=Ftemp(m);
                        COUNT(n)=COUNT(n)+1;
                        
                        
                    end
                end
                k1=k1+1;
            end
            
            for h=2:MAXharmonic
                dBI{h-1} = interp1(Flog{h},dB{h},FT{2});
                POWER=POWER+10.^(dBI{h-1}'/10);    
            end
            Fpower=FT{2};
            
            
            clear ktemp Ftemp COUNT I
            
            % figure,semilogx(...
            %     F1,20*log10(abs(Fpart1(1:end/2)))+CAL,'o:',...
            %     F1log,dB1+CAL,'k*',...
            %     F2,20*log10(abs(Fpart2(1:end/2)))+CAL,'ro',...
            %     F2log,dB2+CAL,'r*',...
            %     F3,20*log10(abs(Fpart3(1:end/2)))+CAL,'bo',...
            %     F3log,dB3+CAL,'b*',...
            %     F4,20*log10(abs(Fpart4(1:end/2)))+CAL,'g',...
            %     F4log,dB4+CAL,'g*','linewidth',1);
            
            % I1=find(and(F1log>=f_start,F1log<=f_end));
            % I2=find(and(F2log>=f_start,F2log<=f_end/2));
            % I3=find(and(F3log>=f_start,F3log<=f_end/3));
            % I4=find(and(F4log>=f_start,F4log<=f_end/4));
            
            INDEXmedian=find(and(FT{1}>=Konstant_level_start,FT{1}<=Konstant_level_end));
            MEDIAN=median(dBT{1}(INDEXmedian));
            DISTLEVEL(LOOP)=MEDIAN;
            AVVIKELSE_MINUS(LOOP)=abs(max(dBT{1}(INDEXmedian))-MEDIAN);
            AVVIKELSE_PLUS(LOOP)=abs(MEDIAN-min(dBT{1}(INDEXmedian)));
            %Beräkna dist i procent
            POW=sqrt(POWER(1:length(Fpower)))./ (10.^(DISTLEVEL(LOOP)/20));
            DIST_per{LOOP}=POW./(1+POW);
            
            PALETT{1}='r:';
            PALETT{2}='r';
            PALETT{3}='b:';
            PALETT{4}='b';
            PALETT{5}='g:';
            PALETT{6}='g';
            PALETT{7}='c:';
            PALETT{8}='c';
            PALETT{9}='m:';
            PALETT{10}='m';        
            PALETT{11}='y:';
            PALETT{12}='y';        
            
            
            figure,
            %Plotta andra ton
            semilogx(Flog{2},dB{2}-DISTLEVEL(LOOP),PALETT{1},FT{2},dBT{2}-DISTLEVEL(LOOP),PALETT{2},'linewidth',2)
            STRI{1}=sprintf('2:nd harmonic');
            STRI{2}=sprintf('Truncated');
            STRI{3}=sprintf('3:rd harmonic');
            STRI{4}=sprintf('Truncated');
            
            
            hold on
            for k=2:MAXharmonic-1
                semilogx(Flog{k+1},dB{k+1}-DISTLEVEL(LOOP),PALETT{2*k-1},FT{k+1},dBT{k+1}-DISTLEVEL(LOOP),PALETT{2*k},'linewidth',2);
                
                if(k>=3)
                    STRI{2*k-1}=sprintf('%d:th harmonic',k+1);
                    STRI{2*k}=sprintf('Truncated',k+1);
                end
                
            end
            %Plotta THD
            semilogx(Fpower,20*log10(POW),'k','linewidth',2)
            STRI{2*k+1}=sprintf('THD');
            
            grid on
            
            STRING=sprintf('Trolldist v. %3.2f Beta - Dist. at constant SPL ~ %3.1f dB| %s',Trollver,DISTLEVEL(LOOP)+CAL,ELEMENT);
            title(STRING)
            xlim([Konstant_level_start Konstant_level_end]);
            ylim([-80 0]);
            xlabel('Frequency [Hz]');
            STRINGY=sprintf('Dist [dB] @ %d cm',MIKAVSTAND);
            ylabel(STRINGY);
            if(TICK==1)
                set(gca,'XTick',XTICK);
                set(gca,'XTickLabel',XTICK);
            end
            
            legend(STRI);
            drawnow
            
            %Spara plot till fil
            STRING=sprintf('%s%s_DIST_at_%ddB.bmp',PATH,ELEMENT,Konstant_level_LEVEL(LOOP));
            disp(['Saving plot to file: ' STRING]);
            print('-dbmp16m',STRING) 
            
            
            %clear I1 I2 I3 I4 F2 F3 F4 Fpart1 Fpart2 Fpart3 Fpart4 i k m
            
        end
        
        %Färgpalett för plottningen
        COLOR{1}='r';
        COLOR{2}='b';
        COLOR{3}='g';
        COLOR{4}='c';
        COLOR{5}='m';
        COLOR{6}='r--';
        COLOR{7}='b--';
        COLOR{8}='g--';
        COLOR{9}='c--';
        COLOR{10}='m--';
        COLOR{11}='r-.';
        COLOR{12}='b-.';
        COLOR{13}='g-.';
        COLOR{14}='c-.';
        COLOR{15}='m-.';
        COLOR{16}='k';
        H=figure;loglog(Fpower,100*DIST_per{1},COLOR{1},'linewidth',2);
        hold on
        k=1;
        STR{k}=sprintf('SPL = %2.1f dB (+%2.1f -%2.1f dB)',DISTLEVEL(k)+CAL,AVVIKELSE_PLUS(k),AVVIKELSE_MINUS(k));
        for k=2:LOOP
            if(k<=length(COLOR))
                kk=k;
            end
            STR{k}=sprintf('SPL = %2.1f dB (+%2.1f -%2.1f dB)',DISTLEVEL(k)+CAL,AVVIKELSE_PLUS(k),AVVIKELSE_MINUS(k));
            loglog(Fpower,100*DIST_per{k},COLOR{kk},'linewidth',2)
        end
        hold off
        grid on
        STRING=sprintf('THD (%d harmonics) in %% at constant SPL at %d cm | %s ',MAXharmonic,MIKAVSTAND,ELEMENT);
        title(STRING)
        legend(STR)
        xlim([Konstant_level_start Konstant_level_end]);
        ylim([0.03 30]);
        STRING=sprintf('Frequency [Hz] | Trolldist ver. %2.2f',Trollver);
        xlabel(STRING);
        STRINGY=sprintf('Dist in %% @ %d cm',MIKAVSTAND);
        ylabel(STRINGY);
        set(gcf,'Position',[100 100 480 480]);
        YTICK=[0.03;0.1;0.3;1;3;10;30];             %Dist i %
        if(TICK==1)
            set(gca,'XTick',XTICK);
            set(gca,'XTickLabel',XTICK);
        end
        set(gca,'YTick',YTICK);
        set(gca,'YTickLabel',YTICK);
        set(gca,'YMinorGrid','off')
        set(gca,'YMinorTick','off')
        
        drawnow
        
        %Spara plot till fil
        STRING=sprintf('%s%s_THD.bmp',PATH,ELEMENT);
        disp(['Saving plot to file: ' STRING]);
        print('-dbmp16m',STRING) 
        
        %     %Tillväxt
        %     figure;
        %     START_AT=7;
        %     for k=START_AT:LOOP-1
        %         if(k<=length(COLOR))
        %             kk=k;
        %         end
        %         if(k==START_AT+1)
        %             hold on
        %         end
        %         STR{k}=sprintf('SPL = %3.1f dB (+%2.1f -%2.1f dB)',DISTLEVEL(k)+CAL,AVVIKELSE_PLUS(k),AVVIKELSE_MINUS(k));
        %         loglog(Fpower,100*DIST_per{k+1}./(100*DIST_per{START_AT}),COLOR{kk},'linewidth',2)
        %     end
        %     hold off
        %     grid on
        %     STRING=sprintf('THD ökning (i ggr) jämfört med mätning vid %3.1f dB (SPL) | %s',DISTLEVEL(START_AT)+CAL,ELEMENT);
        %     title(STRING)
        %     legend(STR(START_AT+1:k+1))
        %     xlim([Konstant_level_start Konstant_level_end]);
        %     xlabel('Frequency [Hz]');
        %     xlim([Konstant_level_start Konstant_level_end]);
        %     ylim([1 10]);
        
        
    end
    STRING=sprintf('THANK YOU for using the freeware Trolldist ver. %2.2f',Trollver);
    disp(STRING);
    disp('PRESS ENTER WHEN READY - MATLAB WILL CLOSE DOWN (MATLAB Problem with the Soundcard)');
    pause
    exit
    
    %Så var allt klart