Mäta dist

[lilltroll]

Hur ser impulssvaren ut?

[Svante]

Svante, har du provat att mäta hemma med trolldist ?

Sakta i backarna nu... En annan har ju precis skaffat matlab på jobbet.

[phon]

Hur ser impulssvaren ut?

Dom har jag ingen ordning på (heller ... ), men så här bra är ett gammalt Philips dome mellanregister:


Och så här ser det samtidigt ut med SweDist:



Samma okalibrerade alltihop, det är mätt lite snabbt när jag passerade jobbet idag. Kul att det visar nåt i alla fall.

Undrar hur starkt det är egentligen? Jag hade nån volt signal och mätte runt 5 cm framför membranet med en halvkass mik och ljudkortets mikingång.

Kalibreringen på Trolldist jobbar i alla fall såg jag. Jag har ändrat tonen till 1kHz. Lite ändrade skalor också, och startfrekv, stoppfrekv, samplingsfrekv plus nåt mer, fick en radda felmeddelanden på slutet. Det här är ju kul

[lilltroll]

Ja, det är nog mycket som kan gå fel i Trolldist, för jag har inte skrivit någon kod som tar hand om / varnar när något är fel. Det får Svante göra om han portar det till BASTA. Personligen tycker jag att MATLAB är bra till prototyper.

Kul att du fick upp kurvor, och att någon utomstående provar! Jag är inte helt med på att du fick det att fungera med bara mikrofoningång.
MATLAB spelar in i stereo i trolldist. Kommer det samma signal i vänster och höger från ljudkortet med mikrofonen inkopplad ?

[lilltroll]

Nu kan man trunkera impulssvaret innan första reflexen, som i MLSSA och andra

Någrar dagar till så har jag en prototyp för ett fullfjädrat mätprog, (dock kanske inte så användarvänligt )

[lilltroll]

Hoppas ni inte blir olyckliga nu ... ännu fler val
Manual kan ni glömma

PS. Man ska trycka 4 ggr med musknapparna för att trunkera fönstret i växande tidsordning. DS

Kod: Markera allt

%Trolldist BETA v 3.94
%Dist & frequency response measuring program using log-swept sine
%Mikael Bohman VT06
%MATLAB Version 7.1.0.246 (R14) Service Pack 3
close all
clear all
clc


fs=44100;   %Samplingsfrekvens [Hz]
T=10;       %Sveptid [s]
f_start=20;    %Startfrekvens [Hz]
f_end=20000;%    %Slutfrekvens [Hz]
CAL=0;          %Kalibrering dB
KONTROLL=1;      % Kontroll av invers
INITTEST=1;     %Initialtest av nivån
FREQ=1000;      %Frekvens vid initaltest [Hz]
Time=2;         %Längden på initialtest [s]

PREDELAY=512;     % Predelay omgång 1 i samples
Prefactor=10;       %Chirpen börjar på en frekvens som är f_start/Precator
Postfactor=2205/2000;     %Chirpen slutar på en frekvens som är f_end/Postfactor
Prewindow=8;        %1/Prefactor fönstras i början
Postwindow=100;      %1/Postwindow fönstras i slutet
Simulate=1;         %Simulera / Mät dist
Konstant_level=0;   %Processen med konstant nivå
Konstant_level_start=28;    %Den frekvens som konstant utnivå börjar på
Konstant_level_end=400;   %Den frekvens som konstant utnivå slutar på
RES=200;                    %Antalet punkters upplösning i frekvensled på logskala
ZOOM=2000;                  %Antalet sampel som visas vid fönstringen av IR

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

if(INITTEST>0)  ;%Initialtest av nivån
    disp('Playing test signal');
    t=0:1/fs:Time;
    SIN=0.99*sin(2*pi*FREQ*t);
    sound(SIN,fs);%,DATA=wavrecord(length(t),fs,2);
    pause(Time+0.3)
    clear SIN t DATA
end

disp('Calculating signal');
t=0:1/fs:T-1/fs;

%y=tdSig-mean(tdSig);
%yinv=tdSigInv-mean(tdSigInv);

y = 0.99*chirp(t,f_start/Prefactor,T,f_end*Postfactor,'logarithmic',-90); %Skapa chirp
n=linspace(-pi,0,ceil(length(y)/Prewindow));
m=linspace(0,pi,ceil(length(y)/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
k=log( Postfactor * f_end ./ f_start *Prefactor)/T;     %Beräknar tidskostanten för inversen
yinv=y(end:-1:1).*exp(-t*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
P=mean(abs(Y.*YINV).^2);                                %Beräkna medeleffekten av Y(f)*YINV(f)
yinv=yinv./sqrt(P);                                     %Korrigera amplituden på inversen

if(KONTROLL>0)
    f=linspace(0,fs/2,length(y));
    figure,semilogx( f, 20*log10(abs(YINV(1:end/2)./sqrt(P)))...
        , f , 20*log10(abs(Y(1:end/2)) )...
        , f , 20*log10( abs( Y(1:end/2).*YINV(1:end/2)./sqrt(P) )));
    xlim([f_start f_end]);
    xlabel('Frequency [Hz]');
    ylabel('Level [dB]');
    legend('Inverted signal i frequency domain','Signal in frequency domain','Product in frequency domain');
   
    drawnow
end
clear Y YINV n m f P

if(Simulate>0)
    %Mät dist
    disp('Recording signal')
    sound(y,fs);REC=wavrecord(length(y),fs,2);
    disp('Processing recording');
    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 y
else
  %Simulera dist
    disp('Simulating dist.')
    [b,a]=butter(3,[0.01 0.1]);
    DATA=0.01*y(:)+(filter(b,a,y(:))).^2;
    REF=y(:);
    clear y b a
end

DATA=[DATA' zeros(1,length(DATA))];
REF=[REF' zeros(1,length(REF))];
out=fftfilt(yinv,DATA);
ref=fftfilt(yinv,REF);

[apa, i]=max(abs(ref)); %Jämför med referenssignalen

clear ref REF DATA apa yinv

L1=fs;
L2=ceil(log(2)/k*fs);
L3=ceil(log(3)/k*fs);
L4=ceil(log(4)/k*fs);
INDEX1=i-PREDELAY:i-PREDELAY+L1-1;
INDEX2=i-PREDELAY-L2:i-PREDELAY-L2-1+round(0.5*L2);
INDEX3=i-PREDELAY-L3:i-PREDELAY-L3-1+round(0.5*(L3-L2));
INDEX4=i-PREDELAY-L4:i-PREDELAY-L4-1+round(0.5*(L4-L3));
part1=out(INDEX1).*tukeywin(length(INDEX1),2*PREDELAY/length(INDEX1))';
part2=out(INDEX2).*tukeywin(length(INDEX2),2*PREDELAY/length(INDEX2))';
part3=out(INDEX3).*tukeywin(length(INDEX3),2*PREDELAY/length(INDEX3))';
part4=out(INDEX4).*tukeywin(length(INDEX4),2*PREDELAY/length(INDEX4))';
Fpart1=fft(part1);
Fpart2=fft(part2);
Fpart3=fft(part3);
Fpart4=fft(part4);

clear INDEX1 INDEX2 INDEX3 INDEX4

figure,
subplot(2,2,1),plot(0:1000/fs:1000*(length(part1)-1)/fs,1000*part1);
title('Impulse response (IR), fundamental')
xlabel('Time [ms]')
subplot(2,2,2),plot(0:1000/fs:1000*(length(part2)-1)/fs,1000*part2);
title('IR 2:nd harmonic')
xlabel('Time [ms]')
subplot(2,2,3),plot(0:1000/fs:1000*(length(part3)-1)/fs,1000*part3);
title('IR, 3:nd harmonic')
xlabel('Time [ms]')
subplot(2,2,4),plot(0:1000/fs:1000*(length(part4)-1)/fs,1000*part4);
title('IR, 4:nd harmonic')
xlabel('Time [ms]')
drawnow

H_1=figure,plot(0:1000/fs:1000*(ZOOM-1)/fs,1000*part1(1:ZOOM));
title('IR-fundamental; Press mousebutton 4 times at the four (start-stop start-stop) raised cosine window locations')
xlabel('Time [ms]')
BUTTON=0;
[Xg,Yg,BUTTON]=ginput(4);    %TRYCK 4 GGR MED MUSEN
N=round(Xg*fs/1000);
TRUNK=part1(N(1):N(4));
WINDOW1=hann(2*(N(2)-N(1)))';
TRUNK(1:N(2)-N(1))=TRUNK(1:N(2)-N(1)).*WINDOW1(1:N(2)-N(1));
WINDOW2=hann(2*(N(4)-N(3)))';
TRUNK(N(3)-N(1):N(4)-N(1))=TRUNK(N(3)-N(1):N(4)-N(1)).*WINDOW2(N(4)-N(3):end);
    %TRUNK(1:N(2)-N(1))=TRUNK(1:N(2)-N(1))*0.5.*(1+cos(0:N(2)-N(1))));
   Ftrunk=fft(TRUNK);
   MAX=max(abs(TRUNK));
figure,plot(0:1000/fs:1000*(length(TRUNK)-1)/fs,TRUNK,...
    0:1000/fs:1000*(length(WINDOW1)/2-1)/fs,MAX.*WINDOW1(1:length(WINDOW1)/2),'r:',...
    Xg(3)-Xg(1):1000/fs:Xg(3)-Xg(1)+1000*(length(WINDOW2)/2)/fs,MAX.*WINDOW2(length(WINDOW2)/2:end),'r:' )
legend('Impulse','Raised cosine windows')
xlabel('Time [ms]')
title('Truncated impulse response')
%clear part1 part2 part3 part4

t=-i/fs:1/fs:(length(out)-i-1)/fs ;
figure, plot(t,10*log10(out.^2));
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);
plot([-L2/fs -L2/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
plot([-L3/fs -L3/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
plot([-L4/fs -L4/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',2);
hold off
legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th harmonic','Location','NorthWest');
drawnow
clear L1 L2 L3 L4 t out


F1=linspace(0,fs/2,length(Fpart1)/2);
F2=linspace(0,fs/4,length(Fpart2)/2);
F3=linspace(0,fs/6,length(Fpart3)/2);
F4=linspace(0,fs/8,length(Fpart4)/2);
FT=linspace(0,fs/2,length(Ftrunk)/2);

[apa, I]=find(F1>1000);
clear apa
CAL=100-20*log10(abs(Fpart1(I(1)))); % Kalibrering

ktemp=log(f_end/f_start);
Ftemp=f_start*exp(ktemp*linspace(0,1,RES));
k1=1;k2=1;k3=1;k4=1;
for m=2:2:length(Ftemp)-1
   
    I=find(and(F1>Ftemp(m-1), F1<Ftemp(m+1)));
    if(isempty(I)~=1)
        dB1(k1)=10*log10(mean(abs(Fpart1(I)).^2));
        F1log(k1)=Ftemp(m);
        k1=k1+1;
    end
    I=find(and(F2>Ftemp(m-1), F2<Ftemp(m+1)));
    if(isempty(I)~=1)
        dB2(k2)=10*log10(mean(abs(Fpart2(I)).^2));
        F2log(k2)=Ftemp(m);
        k2=k2+1;
    end
    I=find(and(F3>Ftemp(m-1), F3<Ftemp(m+1)));
    if(isempty(I)~=1)
        dB3(k3)=10*log10(mean(abs(Fpart3(I)).^2));
        F3log(k3)=Ftemp(m);
        k3=k3+1;
    end
    I=find(and(F4>Ftemp(m-1), F4<Ftemp(m+1)));
    if(isempty(I)~=1)
        dB4(k4)=10*log10(mean(abs(Fpart4(I)).^2));
        F4log(k4)=Ftemp(m);
        k4=k4+1;
    end
end

dBT=20*log10(abs(Ftrunk(1:length(Ftrunk)/2)));
clear ktemp Ftemp k1 k2 k3 k4 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));

figure,semilogx(...
    F1log(I1),dB1(I1)+CAL,'k:',...
    FT,dBT+CAL,'k',...
    F2log(I2),dB2(I2)+CAL,'r:',...
    F3log(I3),dB3(I3)+CAL,'b:',...
    F4log(I4),dB4(I4)+CAL,'g:','linewidth',2);

grid on

xlim([f_start f_end]);
%ylim([10 110]);
xlabel('Frequency [Hz]');
ylabel('SPL [dB] @ 1m ref 20 uPa');
legend('Fundamental','Truncated before first reflex','2:nd harmonic','3:d harmonic','4:th harmonic','Location','NorthWest');
drawnow
title('Trolldist v. 3.94 Beta - Frequency resp. and dist.')
clear I1 I2 I3 I4 F2 F3 F4 Fpart1 Fpart2 Fpart3 Fpart4 i k m

[lilltroll]

PHON, med RES=XXX (som är lite kryptisk för tillfället) kan du ställa upplösningen i frekvensled i senaste Trolldist ver. Den kanske var i grövsta laget ovan.

[phon]

Hjälp

Det är ju dubbelt så många rader nu, mycket att läsa. Manualen försöker man fundera ut on the fly, så att säga ... det är ju halva nöjet. Jodå, kurvor blir det, men jag har inte kommit längre än till distfönstret än ....

Jo, miken ger nog signal till både höger och vänster kanal, ljudet ut kör jag bara ut i en kanal när jag mäter. Det är en liten skräpelektret, men moddad till tretråds source-följare, och en liten skrivbordsförstärkare IC på några w. Lagom stort att leka med mellan arbetsuppgifterna


Ta hand om fel? Det gör ju Matlab åt en, stannar och ger felmeddelanden, fungerar väl bra?




Edit:

så här ser det ut idag, tjusig va?



sedan stannade alltihop, typ såhär:

Kod: Markera allt

H_1 =

     3

??? Error using ==> c:/matlab6p5/toolbox/signal/signal/private/check_order
Order cannot be less than zero.

Error in ==> C:\MATLAB6p5\toolbox\signal\signal\private\gencoswin.m
On line 18  ==> [n,w,trivialwin] = check_order(n);

Error in ==> C:\MATLAB6p5\toolbox\signal\signal\hann.m
On line 19  ==> [w,msg] = gencoswin('hann',varargin{:});

Error in ==> C:\MATLAB6p5\work\Trolldist_394.m
On line 150  ==> WINDOW2=hann(2*(N(4)-N(3)))';

>>

[lilltroll]

Provade att köra koden i MATLAB 6.5

COMSOL har ändrat rutinen för chirp.

En log-chirp i ver 6.5 ÄR INTE samma sak som en log-chirp i ver 7.X

Därför räknar Trolldist FEL under MATLAB 6.5.

Mycket skrämmande

(6.5 verkar beräkna log-chirpen fel!)

En ny Trolldist är på väg som fungerar från MATLAB 6.5 och framåt

[lilltroll]

Nu har trolldist nått nya höjder. Har skrivit om ganska mycket kod, PHON hoppas du inte är rädd för cell i MATLAB

Två mätningar i LEGO-låda på diskanter.
Allt är identiskt förutom bytet av diskant.

SPL skalan är kalibrerad. Exitationsignalen har en RMS amplitud på 2.83 V.

B&K mikrofon samt USB- Creative Extigy används TILLSAMMANS MED BETA VER. 5.7 vilken har en bugg, vilket medför att
kurvan tappar över 14 kHz. Snart kommer byte till buggfriare ver.





VILKEN SKILLNAD runt 5-6 k Hz Haakan hade rätt !

[lilltroll]

Om vi nu jämför de båda diskanterna vid 94 dB ljudtrycksnivå nivå ref 20uPa @ 50 cm.

SCANSPEAK 9300 i LEGO låda. Trolldist v5.8


Peerless-vline_DX25TG05-04.GIF i LEGO låda. Trolldist v5.8

[lilltroll]

Mera mördarkod
>500 rader nu

Kod: Markera allt

%Trolldist BETA v 5.8
%Dist & frequency response measuring program using log-swept sine
%Mikael Bohman VT06
%MATLAB Version 7.1.0.246 (R14) Service Pack 3
%MATLAB Version 6.5.0 (R13)

close all
clear all
clc


fs=48000;       %Samplingsfrekvens [Hz]
T=10;           %Sveptid [s]
AMP=1/1.5;      %Amplitud på utsignalen <1
f_start=100;    %Startfrekvens [Hz]
f_end=20000;%   %Slutfrekvens [Hz]
FREQ=2000;      %Frekvens vid initaltest [Hz]
Time=1;         %Längden på initialtest [s]
CALdB=94;       %SPL nivå vid kalibrering [dB]
CALFRERQ=7000;  %Kalibreringsfrekvens [Hz]

PREDELAY=128;       %Predelay omgång 1 i samples
Prefactor=2;        %Chirpen börjar på en frekvens som är f_start/Precator
Postfactor=0.5*fs/f_end;    %Chirpen slutar på en frekvens som är f_end/Postfactor
Prewindow=8;            %1/Prefactor fönstras i början
Postwindow=100;         %1/Postwindow fönstras i slutet
Konstant_level_start=800;    %Den frekvens som konstant utnivå börjar på
Konstant_level_end=10000;     %Den frekvens som konstant utnivå slutar på
RES=200;                    %Antalet punkters upplösning i frekvensled på logskala
ZOOM=round(fs*0.02);        %Antalet sampel som visas vid fönstringen av IR
MAXharmonic=4;              %Hur många övertoner beräknas totalt (Påverkar THD)

INITTEST=1;     %Initialtest av nivån
KONTROLL=1;      % Kontroll av invers
Konstant_level=1;   %Processen med konstant nivå
Simulate=1;         %Simulera / Mät dist

ELEMENT='SCANSPEAK 9300';

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

if(INITTEST>0)  ;%Initialtest av nivån
    disp('Playing test signal');
    t=0:1/fs:Time;
    SIN=0.99*sin(2*pi*FREQ*t);
    sound(SIN,fs);%,DATA=wavrecord(length(t),fs,2);
    pause(Time+0.3)
    clear SIN t DATA
end

disp('Calculating signal');
t=0:1/fs:T-1/fs;

%y=tdSig-mean(tdSig);
%yinv=tdSigInv-mean(tdSigInv);

%y = 0.99*chirp(t,,T,f_end*Postfactor,'logarithmic',-90); %Skapa chirp
f0=f_start/Prefactor;t1=T;f1=f_end*Postfactor;phi=-90;
instPhi = t1/log(f1/f0)*(f0*(f1/f0).^(t/t1)-f0);
y=0.99*cos(2*pi*(instPhi+phi/360));

clear instPhi

n=linspace(-pi,0,ceil(length(y)/Prewindow));
m=linspace(0,pi,ceil(length(y)/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
time_k=log( Postfactor * f_end ./ f_start *Prefactor)/T;    %Beräknar tidskostanten för inversen
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
P=mean(abs(Y.*YINV).^2);                                %Beräkna medeleffekten av Y(f)*YINV(f)
yinv=yinv./sqrt(P);                                     %Korrigera amplituden på inversen

if(KONTROLL>0)
    f=linspace(0,fs/2,length(y));
    figure,semilogx( f, 20*log10(abs(YINV(1:end/2)./sqrt(P)))...
        , f , 20*log10(abs(Y(1:end/2)) )...
        , f , 20*log10( abs( Y(1:end/2).*YINV(1:end/2)./sqrt(P) )));
    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
end
clear Y YINV n m f P

if(Simulate>0)
    %Mät dist
    disp('Recording signal')
    sound(AMP*y,fs);REC=wavrecord(length(y),fs,2);
    disp('Processing recording');
   
    if(max(max(abs(REC)>0.9)))
        disp('!WARNING!, input signal is overloaded\nProgram halted')
        break
    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
else
    %Simulera dist
    disp('Simulating dist.')
    [b,a]=butter(3,[0.01 0.1]);
    DATA=y(:)+0.01*filter(b,a,y(:)).^2;
    REF=y(:);
    clear b a
end

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)

clear ref REF DATA apa

L{1}=fs;
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
figure,
for n=1:MAXharmonic
    part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
    Fpart{n}=fft(part{n});
    subplot(2,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

clear INDEX
drawnow

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;
N=zeros(4,1);
%N=round(fs/1000*[3 4 8 9])

n=1;
while( ( N(1)<N(2)&& N(2)<N(3) && N(3)<N(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=round(Xg*fs/1000);
    n=n+1;
end

WINDOW1=hann(2*(N(2)-N(1)))';
WINDOW2=hann(2*(N(4)-N(3)))';

for n=1:MAXharmonic
    TRUNK{n}=part{n}(N(1):N(4));
    TRUNK{n}(1:N(2)-N(1))=TRUNK{n}(1:N(2)-N(1)).*WINDOW1(1:N(2)-N(1));
    TRUNK{n}(N(3)-N(1):N(4)-N(1))=TRUNK{n}(N(3)-N(1):N(4)-N(1)).*WINDOW2(N(4)-N(3):end);
    Ftrunk{n}=fft(TRUNK{n});
end
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(3)*1000/fs-N(1)*1000/fs:1000/fs:N(3)*1000/fs-N(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')
%clear part1 part2 part3 part4

t=-i/fs:1/fs:(length(out)-i-1)/fs ;
figure, plot(t,10*log10(out.^2));
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);

plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
hold off
legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th harmonic');
drawnow
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(Ftrunk{n})/2);
    dBT{n}=20*log10(abs(Ftrunk{n}(1:length(Ftrunk{n})/2)));
end


[apa, I]=find(F{1}>CALFRERQ);
clear apa
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(2*RES,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)).^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


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));

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',...
    Fpower,10*log10(POWER(1:length(Fpower))+eps)+CAL,'g:','linewidth',2);

grid on

xlim([f_start f_end]);
ylim([20 100]);
xlabel('Frequency [Hz]');
ylabel('SPL [dB] @ 50 cm ref 20 uPa, 2.83V');
legend('Fundamental','Truncated before first reflex','2:nd harmonic',...
    'Truncated before first reflex','3:d harmonic','Truncated before first reflex','THD+N');
drawnow
title(['Trolldist v. 5.8 Beta - Frequency resp. and dist. | ' ELEMENT])


if(Konstant_level>0)
   
    %Konstant utnivå
   
    [APA,Il]=find(FT{1}<Konstant_level_start); %Hz
    [APA,Ih]=find(FT{1}>Konstant_level_end); %Hz
   

    LEVEL=dBT{1}-dBT{1}(Il(end));
    LEVEL(Il)=LEVEL(Il(end));
    LEVEL(Ih)=LEVEL(Ih(1));
    A=10.^(-LEVEL/20);
   
    %A(Il)=linspace(0,1,length(A(Il)));              %Triangelfönster utanför området
    %A(Ih)=A(Ih(1)).*linspace(1,0,length(A(Ih)));    %Triangelfönster utanför området
    B=fir2(fs,FT{1}./(0.5*fs),A);
   
    yB=fftfilt(B,[y zeros(1,fs)]);
    %clear y
    yB=yB(fs/2:end-fs/2-1);
   
    disp('Recording signal')
    soundsc(yB,fs);REC=wavrecord(length(yB),fs,2);
    disp('Processing recording');
   
    if(max(max(abs(REC)>0.9)))
        disp('!WARNING!, input signal is overloaded\nProgram halted')
        break
    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 y
   
    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)
   
    clear ref REF DATA apa yinv
   
    L{1}=fs;
    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
    figure,
    for n=1:MAXharmonic
        part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
        Fpart{n}=fft(part{n});
        subplot(2,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
   
    clear INDEX
    drawnow
   
    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;
    N=zeros(4,1);
    N=round(fs/1000*[3 4 8 9])
   
    n=1;
    while( ( N(1)<N(2)&& N(2)<N(3) && N(3)<N(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=round(Xg*fs/1000);
        n=n+1;
    end
   
    WINDOW1=hann(2*(N(2)-N(1)))';
    WINDOW2=hann(2*(N(4)-N(3)))';
   
    for n=1:MAXharmonic
        TRUNK{n}=part{n}(N(1):N(4));
        TRUNK{n}(1:N(2)-N(1))=TRUNK{n}(1:N(2)-N(1)).*WINDOW1(1:N(2)-N(1));
        TRUNK{n}(N(3)-N(1):N(4)-N(1))=TRUNK{n}(N(3)-N(1):N(4)-N(1)).*WINDOW2(N(4)-N(3):end);
        Ftrunk{n}=fft(TRUNK{n});
    end
    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(3)*1000/fs-N(1)*1000/fs:1000/fs:N(3)*1000/fs-N(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')
    %clear part1 part2 part3 part4
   
    t=-i/fs:1/fs:(length(out)-i-1)/fs ;
    figure, plot(t,10*log10(out.^2));
    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);
   
    plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
    plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
    plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
    hold off
    legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th harmonic');
    drawnow
    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(Ftrunk{n})/2);
        dBT{n}=20*log10(abs(Ftrunk{n}(1:length(Ftrunk{n})/2)));
    end
   
   
    [apa, I]=find(F{1}>CALFRERQ);
    clear apa
    %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(2*RES,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)).^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
   
   
    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));
   
    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',...
        Fpower,10*log10(POWER(1:length(Fpower))+eps)+CAL,'g:','linewidth',2);
   
    grid on
   
    xlim([Konstant_level_start Konstant_level_end]);
    ylim([20 100]);
    xlabel('Frequency [Hz]');
    ylabel('SPL [dB] @ 50 cm ref 20 uPa, 2.83V');
    legend('Fundamental','Truncated before first reflex','2:nd harmonic',...
        'Truncated before first reflex','3:d harmonic','Truncated before first reflex','THD+N');
    drawnow
    title(['Trolldist v. 5.8 Beta - Dist. at constant SPL ~ 94 dB| ' ELEMENT])
    clear I1 I2 I3 I4 F2 F3 F4 Fpart1 Fpart2 Fpart3 Fpart4 i k m
   
end

[phon]

Oj, såg inte nya koden Har precis testat 5.7 och den gick i stort sett bra. Har Matlab 6.5.0 R13 precis som det står i koden, det finns tydligen bättre.

Har inget att svepa på och med där jag sitter, bara en 90 dB väsnandes PC med hörselskador.


Nu skall jag prova 5.8


edit:

å den gick också, fast samma felkod.

Är det min log10.m som är konstig?

[lilltroll]

Oj, såg inte nya koden Har precis testat 5.7 och den gick i stort sett bra. Har Matlab 6.5.0 R13 precis som det står i koden, det finns tydligen bättre.

Har inget att svepa på och med där jag sitter, bara en 90 dB väsnandes PC med hörselskador.


Nu skall jag prova 5.8


edit:

å den gick också, fast samma felkod.

Är det min log10.m som är konstig?

Log(0) = -Inf, det är det som är varningen. Inget att bry sig om.

[lilltroll]

Nu med THD i %

Kod: Markera allt

%Trolldist BETA v 5.9
%Dist & frequency response measuring program using log-swept sine
%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=5.9;   %Ver.
fs=48000;       %Samplingsfrekvens [Hz]
T=4;            %Sveptid [s]
AMP=1;          %Amplitud på utsignalen vid frekvensgångsmätning AMP<1
AMP2=1;         %Amplitud på utsignalen vid distmätning AMP<1
f_start=100;    %Startfrekvens [Hz]
f_end=20000;%   %Slutfrekvens [Hz]
FREQ=2000;      %Frekvens vid initaltest [Hz]
Time=1;         %Längden på initialtest [s]
CALdB=94.3;       %SPL nivå vid kalibrering [dB]
CALFRERQ=7078;  %Kalibreringsfrekvens [Hz]

WINDOW=[eps 1 5.5 6.5];   %Fönstringsval [ms]

PREDELAY=0;       %Predelay omgång 1 i samples
Prefactor=2;        %Chirpen börjar på en frekvens som är f_start/Precator
Postfactor=0.5*fs/f_end;    %Chirpen slutar på en frekvens som är f_end/Postfactor
Prewindow=8;            %1/Prefactor fönstras i början
Postwindow=200;         %1/Postwindow fönstras i slutet
Konstant_level_start=800;    %Den frekvens som konstant utnivå börjar på
Konstant_level_end=10000;     %Den frekvens som konstant utnivå slutar på
Konstant_level_LEVEL=95;       %Den önskade nivån vid dist mätning
RES=200;                    %Antalet punkters upplösning i frekvensled på logskala
ZOOM=round(fs*0.02);        %Antalet sampel som visas vid fönstringen av IR
MAXharmonic=4;              %Hur många övertoner beräknas totalt (Påverkar THD)

USE_WINDOW=0;   %Använd manuell musfönstring / manuellt förval
INITTEST=1;     %Initialtest av nivån
KONTROLL=1;      % Kontroll av invers
Konstant_level=1;   %Processen med konstant nivå
Simulate=1;         %Simulera / Mät dist

ELEMENT='SCANSPEAK 9300';

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

if(INITTEST>0)  ;%Initialtest av nivån
    disp('Playing test signal');
    t=0:1/fs:Time;
    SIN=0.99*AMP*sin(2*pi*FREQ*t);
    sound(SIN,fs);%,DATA=wavrecord(length(t),fs,2);
    pause(Time+0.3)
    clear SIN t DATA
end

disp('Calculating signal');
t=0:1/fs:T-1/fs;

%y=tdSig-mean(tdSig);
%yinv=tdSigInv-mean(tdSigInv);

%y = 0.99*chirp(t,,T,f_end*Postfactor,'logarithmic',-90); %Skapa chirp
f0=f_start/Prefactor;t1=T;f1=f_end*Postfactor;phi=-90;
instPhi = t1/log(f1/f0)*(f0*(f1/f0).^(t/t1)-f0);
y=0.99*cos(2*pi*(instPhi+phi/360));

clear instPhi

n=linspace(-pi,0,ceil(length(y)/Prewindow));
m=linspace(0,pi,ceil(length(y)/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
time_k=log( Postfactor * f_end ./ f_start *Prefactor)/T;    %Beräknar tidskostanten för inversen
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
P=mean(abs(Y.*YINV).^2);                                %Beräkna medeleffekten av Y(f)*YINV(f)
yinv=yinv./sqrt(P);                                     %Korrigera amplituden på inversen

if(KONTROLL>0)
    f=linspace(0,fs/2,length(y));
    figure,semilogx( f, 20*log10(abs(YINV(1:end/2)./sqrt(P)))...
        , f , 20*log10(abs(Y(1:end/2)) )...
        , f , 20*log10( abs( Y(1:end/2).*YINV(1:end/2)./sqrt(P) )));
    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
end
clear Y YINV f P

if(Simulate>0)
    %Mät dist
    disp('Recording signal')
    sound(AMP*y,fs);REC=wavrecord(length(y),fs,2);
    disp('Processing recording');
   
    if(max(max(abs(REC)>0.9)))
        disp('!WARNING!, input signal is overloaded\nProgram halted')
        break
    end
   
    DATA=REC(:,1);
    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;
    REF=y(:);
    clear b a
end

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)

clear ref REF DATA apa

L{1}=fs;
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
figure,
for n=1:MAXharmonic
    part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
    Fpart{n}=fft(part{n});
    subplot(2,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

clear INDEX
drawnow

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

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);
    Ftrunk{n}=fft(TRUNK{n});
end
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')
%clear part1 part2 part3 part4

t=-i/fs:1/fs:(length(out)-i-1)/fs ;
figure, plot(t,10*log10(out.^2));
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);

plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
hold off
legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th harmonic');
drawnow
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(Ftrunk{n})/2);
    dBT{n}=20*log10(abs(Ftrunk{n}(1:length(Ftrunk{n})/2)));
end


[apa, I]=find(F{1}>CALFRERQ);
clear apa
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(2*RES,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)).^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


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));

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',...
    Fpower,10*log10(POWER(1:length(Fpower))+eps)+CAL,'g:','linewidth',2);

grid on

xlim([f_start f_end]);
ylim([20 100]);
xlabel('Frequency [Hz]');
ylabel('SPL [dB] @ 50 cm ref 20 uPa, 2.83V');
legend('Fundamental','Truncated before first reflection','2:nd harmonic',...
    'Truncated before first reflection','3:rd harmonic','Truncated before first reflection','THD+N');
STRING=sprintf('Trolldist v. %3.1f Beta - Dist. and Frequency resp. | %s',Trollver,ELEMENT);
title(STRING)
set(gcf,'Position',[100 100 480 480]);
drawnow
IR=TRUNK;
time=0:1/fs:(length(IR{1})-1)/fs;
save(ELEMENT,'IR','fs','time')

if(Konstant_level>0)
   
    %Konstant utnivå
   
    [APA,Il]=find(FT{1}<Konstant_level_start); %Hz
    [APA,Ih]=find(FT{1}>Konstant_level_end); %Hz
   

    LEVEL=dBT{1}-dBT{1}(Il(end));
    LEVEL(Il)=LEVEL(Il(end));
    LEVEL(Ih)=LEVEL(Ih(1));
    A=10.^(-LEVEL/20);
   
    %A(Il)=linspace(0,1,length(A(Il)));              %Triangelfönster utanför området
    %A(Ih)=A(Ih(1)).*linspace(1,0,length(A(Ih)));    %Triangelfönster utanför området
    B=fir2(fs,FT{1}./(0.5*fs),A);
   
    yB=fftfilt(B,[y zeros(1,fs)]);
    %clear y
    yB=AMP2*0.99*yB(fs/2:end-fs/2-1)./max(abs(yB(fs/2:end-fs/2-1)));
   
    disp('Recording signal')
    sound(yB,fs);REC=wavrecord(length(yB),fs,2);
    disp('Processing recording');
   
    if(max(max(abs(REC)>0.9)))
        disp('!WARNING!, input signal is overloaded\nProgram halted')
        break
    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 y
   
    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)
   
    clear ref REF DATA apa yinv
   
    L{1}=fs;
    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
    figure,
    for n=1:MAXharmonic
        part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
        Fpart{n}=fft(part{n});
        subplot(2,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
   
    clear INDEX
    drawnow
   
    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;
   
    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);
        Ftrunk{n}=fft(TRUNK{n});
    end
    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')
    %clear part1 part2 part3 part4
   
    t=-i/fs:1/fs:(length(out)-i-1)/fs ;
    figure, plot(t,10*log10(out.^2));
    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);
   
    plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
    plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
    plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
    hold off
    legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th harmonic');
    drawnow
    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(Ftrunk{n})/2);
        dBT{n}=20*log10(abs(Ftrunk{n}(1:length(Ftrunk{n})/2)));
    end
   
   
    [apa, I]=find(F{1}>CALFRERQ);
    clear apa
    %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(2*RES,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)).^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
   
   
    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));
    DISTLEVEL=median(dBT{1}(INDEXmedian));
    POW=sqrt(POWER(1:length(Fpower)))./ (10.^(DISTLEVEL/20));
    DIST_per=POW./(1+POW);
   

   
    figure,
    subplot(2,1,1),semilogx(...
        Flog{2},dB{2}-DISTLEVEL,'r:',...
        FT{2},dBT{2}-DISTLEVEL,'r',...
        Flog{3},dB{3}-DISTLEVEL,'b:',...
        FT{3},dBT{3}-DISTLEVEL,'b',...
        Fpower,10*log10(POWER(1:length(Fpower))+eps)-DISTLEVEL,'g:','linewidth',2);
   
    grid on
   
    STRING=sprintf('Trolldist v. %3.1f Beta - Dist. at constant SPL ~ %3.1f dB| %s',Trollver,DISTLEVEL+CAL,ELEMENT)
    title(STRING)
    xlim([Konstant_level_start Konstant_level_end]);
    ylim([-80 0]);
    xlabel('Frequency [Hz]');
    ylabel('Dist [dB] @ 50 cm');
    legend('2:nd harmonic',...
        'Truncated before first reflex','3:d harmonic','Truncated before first reflex','THD+N [dB]');
   
    subplot(2,1,2),loglog(Fpower,100*DIST_per,'r','linewidth',2)
    grid on
    STRING=sprintf('THD+N in %% at constant SPL~%3.1f dB| %s',DISTLEVEL+CAL,ELEMENT)
    title(STRING)
    xlim([Konstant_level_start Konstant_level_end]);
    ylim([0.1 100]);
    xlabel('Frequency [Hz]');
    ylabel('Dist in % @ 50 cm');
    set(gcf,'Position',[100 100 480 640]);
    drawnow
   
    %clear I1 I2 I3 I4 F2 F3 F4 Fpart1 Fpart2 Fpart3 Fpart4 i k m
   
end

[lilltroll]

Nu Ver 6.0, med en fungerande kalibrering samt en helautomatisk distmätdel

[Kaffekoppen]

[lilltroll]

Fungerande kalibreringsdel
Fungerande distmätningsdel med konstant utnivå vid korta svep
(Baserat på minimum-fas invers EQ m h a Hilbert transform )
Trollets favorit

Beskrivande text i en del av koden

Sparar plottar till fil

Mycket färre buggar

Ingen manual

Kod: Markera allt

%Trolldist BETA v 6.0
%Dist & frequency response measuring program using log-swept sine
%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=6.0;  %Ver.
fs=48000;       %Samplingsfrekvens [Hz]
T=3;            %Sveptid [s]
AMP=10^(-14/20);          %Amplitud på utsignalen vid frekvensgångsmätning AMP<1

%INITIALT TEST MED SINUS
FREQ=1000;      %Frekvens vid initaltest [Hz]
Time=10;        %Längden på initialtest [s]

%KALIBRERINGSSEKVENS
CALdB=86;         %SPL nivå vid kalibrering [dB]
OFFSET_dB=-18.45;   %Offsetvärdet genererat av kalibreringssekvens, Uppdatera denna efter kalibrering!
CALFRERQ=1000;      %Kalibreringsfrekvens [Hz]

%SVEPT SINUS

f_start=100;    %Startfrekvens för svepet [Hz]
f_end=20000;%   %Slutfrekvens för svepet [Hz]
PREDELAY=0;     %Predelay i sample (DEFAULT=0)
Prefactor=2;    %Chirpen börjar på en frekvens som är f_start/Precator
Postfactor=0.5*fs/f_end;    %Chirpen slutar på en frekvens som är f_end/Postfactor
Prewindow=8;                %1/Prefactor fönstras i början
Postwindow=200;             %1/Postwindow fönstras i slutet

%DIST. MÄTNING VID KONSTANT UTNIVÅ

Konstant_level_start=1000;    %Den frekvens som konstant utnivå börjar på
Konstant_level_end=14000;     %Den frekvens som konstant utnivå slutar på
Konstant_level_LEVEL=94;               %Den önskade ljudtrycksnivån vid dist mätning

%DISPLAY

RES=200;                    %Antalet punkters upplösning i frekvensled på logskala
ZOOM=round(fs*0.02);        %Antalet sampel som visas vid fönstringen av IR

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

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

%VAL [0/1]
USE_WINDOW=1;   %Använd ej mus-baserad fönstring = 0 / Använd mus-baserad fönstring = 1
INITTEST=0;     %Initialtest av nivån avslagen = 0 / Initialtest av nivån påslagen=1
KONTROLL=0;     %Grafisk kontroll av inveren = 0 / Grafisk kontroll avslagen = 1
Konstant_level=1;   %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

%Högtalarelement
ELEMENT='53-NDT-811435';

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

%Initialtest av nivån
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

%Beräkna fönster
n=linspace(-pi,0,ceil(length(y)/Prewindow));
m=linspace(0,pi,ceil(length(y)/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
time_k=log( Postfactor * f_end ./ f_start *Prefactor)/T;    %Beräknar tidskostanten för inversen
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
P=mean(abs(Y.*YINV).^2);                                %Beräkna medeleffekten av Y(f)*YINV(f)
yinv=yinv./sqrt(P);                                     %Korrigera amplituden på inversen

%Grafisk kontroll av inversen
if(KONTROLL>0)
    f=linspace(0,fs/2,length(y));
    figure,semilogx( f, 20*log10(abs(YINV(1:end/2)./sqrt(P)))...
        , f , 20*log10(abs(Y(1:end/2)) )...
        , f , 20*log10( abs( Y(1:end/2).*YINV(1:end/2)./sqrt(P) )));
    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
end
clear Y YINV f P

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');
   
    if(max(max(abs(REC)>0.9)))
        disp('!WARNING!, input signal is overloaded\nProgram halted')
        break
    end
   
    DATA=REC(:,1);
    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, i]=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)-i;
if(L{1}>fs)
    L{1}=fs;
end

%Övertonerna
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

%Visa alla estimerade impulsvar

figure,
for n=1:MAXharmonic
    part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
    Fpart{n}=fft(part{n});
    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

drawnow

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

n=1;

%Välj fönstringspunkter m h a musen
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

%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
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_IMP.bmp',ELEMENT);
disp(['Saving plot to file' STRING]);
print('-dbmp16m',STRING)

%Plotta en Energy Time Curve för hela resultatet, samt markera var
%deltonerna 2-4 kan hittas
t=-i/fs:1/fs:(length(out)-i-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);

plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
hold off
legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th harmonic');
drawnow
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

xlim([f_start f_end]);
ylim([20 100]);
xlabel('Frequency [Hz]');
ylabel('SPL [dB] @ 50 cm ref 20 uPa, 2.83V');
legend('Fundamental','Truncated before first reflection','2:nd harmonic',...
    'Truncated before first reflection','3:rd harmonic','Truncated before first reflection'),0;
STRING=sprintf('Trolldist v. %3.1f Beta - Dist. and Frequency resp. | %s',Trollver,ELEMENT);
title(STRING)
set(gcf,'Position',[100 100 480 480]);
drawnow

%Spara plot till fil
STRING=sprintf('%s.bmp',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;
save(ELEMENT,'IR','fs','time')

%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)]))))));
   
    %Beräkna det minimum-fas EQ ade chirpen
    yB=fftfilt(eit,y);
    clear y
   
    %Korrigera till önskad nivå
    AMPdB=Konstant_level_LEVEL-20*log10(abs(Ftrunk{1}(Il(end))))-CAL;
    yB=10.^(AMPdB/20)*0.99*AMP*yB;
   
    %Kontrollera att signalen inte klipper
    MaxA=max(abs(yB));
    if(MaxA>0.99)
        disp('Utvolymen på förstärkaren är inte tillräckligt hög för önskad utnivå')
        disp(MaxA)
        return
    end
   
    %Spela in signalen
    disp('Recording signal')
    pause(0.1);
    wavplay(yB(:)*[1 1],fs,'async');REC=wavrecord(length(yB),fs,2);
    disp('Processing recording');
   
    if(max(max(abs(REC)>0.9)))
        disp('!WARNING!, input signal is overloaded\nProgram halted')
        break
    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 y
   
    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)
   
    clear ref REF DATA apa yinv
   
    %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
    figure,
    for n=1:MAXharmonic
        part{n}=out(INDEX{n}).*tukeywin (length(INDEX{n}),2*PREDELAY/length(INDEX{n}))';
        Fpart{n}=fft(part{n});
        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
   
    clear INDEX
    drawnow
   
    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;
   
    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);
        Ftrunk{n}=fft(TRUNK{n});
    end
    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')
    %clear part1 part2 part3 part4
   
    t=-i/fs:1/fs:(length(out)-i-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);
   
    plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
    plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
    plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
    hold off
    legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th harmonic');
    drawnow
    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(Ftrunk{n})/2);
        dBT{n}=20*log10(abs(Ftrunk{n}(1:length(Ftrunk{n})/2)));
    end
   
   
    [apa, I]=find(F{1}>CALFRERQ);
    clear apa
    %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(2*RES,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)).^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
   
   
    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));
    DISTLEVEL=median(dBT{1}(INDEXmedian));
    %Beräkna dist i procent
    POW=sqrt(POWER(1:length(Fpower)))./ (10.^(DISTLEVEL/20));
    DIST_per=POW./(1+POW);

    figure,
    subplot(2,1,1),semilogx(...
        Flog{2},dB{2}-DISTLEVEL,'r:',...
        FT{2},dBT{2}-DISTLEVEL,'r',...
        Flog{3},dB{3}-DISTLEVEL,'b:',...
        FT{3},dBT{3}-DISTLEVEL,'b','linewidth',2);
%        Fpower(INDEX_FULHACK),10*log10(POWER(1:length(Fpower(INDEX_FULHACK)))+eps)-DISTLEVEL,
   
    grid on
   
    STRING=sprintf('Trolldist v. %3.1f Beta - Dist. at constant SPL ~ %3.1f dB| %s',Trollver,DISTLEVEL+CAL,ELEMENT);
    title(STRING)
    xlim([Konstant_level_start Konstant_level_end]);
    ylim([-80 0]);
    xlabel('Frequency [Hz]');
    ylabel('Dist [dB] @ 50 cm');
    legend('2:nd harmonic','Truncated before first reflex','3:d harmonic','Truncated before first reflex',0);
   
 
    subplot(2,1,2),loglog(Fpower,100*DIST_per,'r','linewidth',2)
    grid on
    STRING=sprintf('THD (%d harmonics) in %% at constant SPL~%3.1f dB| %s',MAXharmonic,DISTLEVEL+CAL,ELEMENT);
    title(STRING)
    xlim([Konstant_level_start Konstant_level_end]);
    ylim([0.1 100]);
    xlabel('Frequency [Hz]');
    ylabel('Dist in % @ 50 cm');
    set(gcf,'Position',[100 100 480 640]);
    drawnow
   
    %Spara plot till fil
    STRING=sprintf('%s_DIST.bmp',ELEMENT);
    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
disp('THANK YOU for using the freeware Trolldist')
disp('PRESS ENTER WHEN READY - MATLAB WILL CLOSE DOWN (MATLAB Problem with the Soundcard)')
pause
exit

%Så var allt klart


[phon]

[lilltroll]

Trolldist 7.0
Jaha, 100 sköna rader längre.

Mäter dist vid flera nivåer om man vill.
Mindre buggar
Mer användarvänlig (Varnar på flertalet saker om något går fel)
Beräknar själv "magiska konstanter" som bara Lilltroll själv kunde ställa i ver <= 6.1
Lite mera hjälptext i koden.
Exporterar koef. för FIR-baserad digital minimumfas EQ - Mycket användbart om man är DSP hacker och vill bygga digitala delningsfilter eller digital EQ
Kanske kommer en superavancerad ver. som kan hitta optimal IIR-baserad EQ
(Jag vet hur man gör Metod 1, m h a Prony's metod 2, minsta kvadratanpassning. Men här tror jag Svante börjar få huvudvärk )
Den kommer smygande redan vid Hilbert transform


Known bugs:
Sparar inte filerna på ett vettigt sätt än. Kommer i ver 7.1

Kod: Markera allt

%Trolldist BETA v 7.0
%Dist & frequency response measuring program using log-swept sine
%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.0;  %Ver.
fs=48000;       %Samplingsfrekvens [Hz]
T=5;            %Sveptid [s] Lämpligt val är inom området 0.5-10 s. T>2*T60 dvs sveptiden bör vara större än 2*efterklangstiden för rummet som man mäter i
AMP=10^(-14/20);          %Amplitud på utsignalen vid frekvensgångsmätning AMP<1. 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=1000;      %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=94.5;         %SPL nivå vid kalibrering [dB], skriv in vad SPL-metern visade här!
OFFSET_dB=-16.2;   %Offsetvärdet genererat av kalibreringssekvens, Uppdatera denna efter kalibrering!
CALFRERQ=1000;      %Kalibreringsfrekvens [Hz]

%SVEPT SINUS

f_start=100;    %Startfrekvens för svepet [Hz]
f_end=20000;%   %Slutfrekvens för svepet [Hz]
PREDELAY=0;     %Predelay i sample (DEFAULT=0)
Prefactor=2;    %Chirpen börjar på en frekvens som är f_start/Precator
Postfactor=1.2;    %Chirpen slutar på en frekvens som är f_end/Postfactor
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=[75 80 85 90 95 100];     %Den önskade ljudtrycksnivånerna vid dist mätning

%DISPLAY

RES=200;                    %Antalet punkters upplösning i frekvensled på logskala
ZOOM=round(fs*0.02);        %Antalet sampel som visas vid fönstringen av IR

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

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

%VAL [0/1]
USE_WINDOW=1;   %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

%Högtalarelement
ELEMENT='CA18RLY'; %Strängen måst4e fungera som filnamn också, inga / här!

disp('Trolldist')
%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('Mätfelet från inversen kommer att ligga inom +- %2.3f dB\nÖnskas ett mindre fel - svep långsammare\n',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');
   
    if(max(max(abs(REC)>0.95)))
        disp('!WARNING!, input signal is overloaded\nProgram halted')
        break
    end
   
    DATA=REC(:,1);
    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, i]=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)-i;
if(L{1}>fs)
    L{1}=fs;
end

%Övertonerna
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

%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


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

n=1;

%Välj fönstringspunkter m h a musen
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

%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_IMP.bmp',ELEMENT);
    disp(['Saving plot to file' STRING]);
    print('-dbmp16m',STRING)
   
    %Plotta en Energy Time Curve för hela resultatet, samt markera var
    %deltonerna 2-4 kan hittas
    t=-i/fs:1/fs:(length(out)-i-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);
   
    plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
    plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
    plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
    hold off
    legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th 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

xlim([f_start f_end]);
ylim([20 100]);
xlabel('Frequency [Hz]');
ylabel('SPL [dB] @ 50 cm ref 20 uPa, 2.83V');
legend('Fundamental','Truncated before first reflection','2:nd harmonic',...
    'Truncated before first reflection','3:rd harmonic','Truncated before first reflection'),0;
STRING=sprintf('Trolldist v. %3.1f Beta - Dist. and Frequency resp. | %s',Trollver,ELEMENT);
title(STRING)
set(gcf,'Position',[100 100 480 480]);
drawnow

%Spara plot till fil
STRING=sprintf('%s.bmp',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')

%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 y 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('Utvolymen på förstärkaren är inte tillräckligt hög för önskad utnivå')
            disp(MaxA)
            return
        end
       
        %Spela in signalen
        disp('Recording signal')
        pause(0.1);
        wavplay(PLAY(:)*[1 1],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 y
       
        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')
            %clear part1 part2 part3 part4
           
            t=-i/fs:1/fs:(length(out)-i-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);
           
            plot([-L{2}/fs -L{2}/fs],[YLIM(2)-100 YLIM(2)],'r','linewidth',2);
            plot([-L{3}/fs -L{3}/fs],[YLIM(2)-100 YLIM(2)],'g','linewidth',2);
            plot([-L{4}/fs -L{4}/fs],[YLIM(2)-100 YLIM(2)],'m','linewidth',1);
            hold off
            legend('Convoluted recording','Fundamental','2:nd harmonic','3:d harmonic','4:th 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
        %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(2*RES,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)).^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
       
       
        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_PLUS(LOOP)=abs(max(dBT{1}(INDEXmedian))-MEDIAN);
        AVVIKELSE_MINUS(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);
       
        figure,
        semilogx(...
            Flog{2},dB{2}-DISTLEVEL(LOOP),'r:',...
            FT{2},dBT{2}-DISTLEVEL(LOOP),'r',...
            Flog{3},dB{3}-DISTLEVEL(LOOP),'b:',...
            FT{3},dBT{3}-DISTLEVEL(LOOP),'b','linewidth',2);
        %        Fpower{LOOP}(INDEX_FULHACK),10*log10(POWER(1:length(Fpower{LOOP}(INDEX_FULHACK)))+eps)-DISTLEVEL,
       
        grid on
       
        STRING=sprintf('Trolldist v. %3.1f 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]');
        ylabel('Dist [dB] @ 50 cm');
        legend('2:nd harmonic','Truncated before first reflex','3:d harmonic','Truncated before first reflex',0);
        drawnow
       
        %Spara plot till fil
        STRING=sprintf('%s_DIST.bmp',ELEMENT);
        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--';
   
   
    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 | %s',MAXharmonic,ELEMENT);
    title(STRING)
    legend(STR)
    xlim([Konstant_level_start Konstant_level_end]);
    ylim([0.05 50]);
    xlabel('Frequency [Hz]');
    ylabel('Dist in % @ 50 cm');
    set(gcf,'Position',[100 100 480 480]);
    drawnow
   
end
disp('THANK YOU for using the freeware Trolldist')
disp('PRESS ENTER WHEN READY - MATLAB WILL CLOSE DOWN (MATLAB Problem with the Soundcard)')
pause
exit

%Så var allt klart

[lilltroll]

Tja

Nu ska vi se så här en lördagsnatt

Alla tidsinvarianta ickelinjära system kan modelleras som en oändlig summa av multidimensionella (faltnings)integraler (påstog Vito Volterra runt 1880)



Den tidsdiskreta motsvarigheten blir:



Trolldist hittar just h1....hn påstår Lilltroll

Det ska således vara möjligt att prediktera alla utsignaler från en godtycklig insignal, baserat på endast en uppsättning h1 .... hn

Men jag förstår inte hur lavindist. kommer in här Gör ni

Ljudåtergivning är nog inte riktigt så enkelt, men borde det inte gå att beskriva det exakt som en oändligt antal faltningar av Volterra kärnor?
Eller, som ett multidimensionellt ickelinjärt system (där varje uppsättning h1....hn tillsammans spänner upp ett Hilbertrum) vilket approx. kan reduceras till en begränsad mängd Volterra kärnor med en begränsat antal h omm systemet är svagt ickelinjärt, (flera dist processer verkar samtidigt, och de måste separeras till varje process för att (enklare) kunna identifiera systemet) ... eller har jag fel där ... för det måste gå att falta flera Volterra kärnor med varandra.
Någon som kan det här, Morello hjälp !

Ett mer konkret exempel: Man kopplar en förstärkare till en högtalare, som genererar ljud i luften. Vet man hur varje del distar så kan man beräkna den totala disten.


Men nu till den svårare frågan:
Hur hittar man en stabil invers till en Volterra kärna?
Dvs hur beräknar jag en stabil ickelinjär utjämnare?

[Svante]

Ujuj, Lilltroll. Tungt...

Inte för att jag är riktigt med i svängarna, men om någon farbror på 1800-talet sa att olinjära system kan beskrivas så, så är det nog riktigt. 1800-talssnubbarna var ju rätt fiffiga.

Men vad du undrar är om det går att göra en invers till ett ickelinjärt system...? Utan att tänka som integralerna du visar, så kan man ju förstå det orimliga i att försöka hitta en inverspryl som försöker återskapa en signal ut från ett klippande slutsteg. Vore slutsteget mjukklippande kanske det skulle gå, nästan, men uj så svårt det blir om man driver det till hårdklippning. Och att hitta en insignal till ett klippande slutsteg, som gör att slutsteget inte klipper, ja det blir snäppet värre.

Är det det du vill göra? Eller vill du framkoppla bort disten i högtalaren, kanske?

[lilltroll]

Ujuj, Lilltroll. Tungt...

Inte för att jag är riktigt med i svängarna, men om någon farbror på 1800-talet sa att olinjära system kan beskrivas så, så är det nog riktigt. 1800-talssnubbarna var ju rätt fiffiga.

Men vad du undrar är om det går att göra en invers till ett ickelinjärt system...? Utan att tänka som integralerna du visar, så kan man ju förstå det orimliga i att försöka hitta en inverspryl som försöker återskapa en signal ut från ett klippande slutsteg. Vore slutsteget mjukklippande kanske det skulle gå, nästan, men uj så svårt det blir om man driver det till hårdklippning. Och att hitta en insignal till ett klippande slutsteg, som gör att slutsteget inte klipper, ja det blir snäppet värre.

Är det det du vill göra? Eller vill du framkoppla bort disten i högtalaren, kanske?

Nej, ett klippande system blir svårt att kompensera. Men då är systemet inte svagt ickelinjärt. Dvs h1-h3 >> h4...hn är inte uppfyllt.

Betänk en enkel olinjäritet. y=x^2

Den kan återskapas med x=sqrt(y), vilken kan serieutvecklas till en lång serie. Men det är just det, det krävs en oändlig serie med övertoner för att kompensera för den enklaste olijäritet genom framkoppling.

Jag kommer bara flytta problemen. Kompenserar jag för 2:a och 3:e ton, så kommer jag få mer dist för de övertoner som jag inte bryr mig om.

Det är högtalarna jag vill kompensera ! Inte förstärkarklippning

[Svante]

Nja, nej, det var kanske mest sett som en illustration till att det i det generella fallet inte alltid går att hitta en invers. Det betyder att man måste specialbehandla varje fall och försäkra sig om att man inte har en klippliknande situation.

[phon]

Det ska således vara möjligt att prediktera alla utsignaler från en godtycklig insignal, baserat på endast en uppsättning h1 .... hn

Men jag förstår inte hur lavindist. kommer in här Gör ni

Prediktiva stokastiska variabler

Det är bara att extrahera ut slumpen ur det hela och mata tillbaka slumpens invers ..... Random-transform?

Alternativt får du väl lägga till en egen slumpgenerator så får du full koll på slumpen.
Det är alltid bra att veta exakt när slumpvisa saker skall inträffa.


EDIT: Är det hysteres på lavineffekten?

[lilltroll]

Trolldist 7.3

860 rader, börjar bli svettigt.

Fixar många små saker. Bättre plottar. Beräknar själv de flesta hemliga konstaterna nu. Bestäm bara sveptid T f_start samt f_end.

Finns fortfarande många förbättringsmöjligheter

Kod: Markera allt

%Trolldist BETA v 7.3
%Dist & frequency response measuring program using log-swept sine
%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.3;  %Ver.

T=5;            %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=-20;       %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=100;         %SPL nivå vid kalibrering [dB], skriv in vad SPL-metern visade här!
OFFSET_dB=-17;   %Offsetvärdet genererat av kalibreringssekvens, Uppdatera denna efter kalibrering!
CALFRERQ=63;      %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=64;     %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=7000;     %Den frekvens som konstant ljudtrycksnivå slutar på
Konstant_level_LEVEL=[80 90];     %Den önskade ljudtrycksnivånerna vid dist mätning

%DISPLAY

RES=150;                    %Antalet punkters upplösning i frekvensled på logskala
ZOOM_ms=30;               %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=[25;30;40;50;60;70;80;90;100;200];    %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.

%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 200 300];   %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=1;   %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!

%Högtalarelement
ELEMENT='XLS'; %Strängen måste fungera som filnamn också, inga / _ ' ' här!
CONDITION='Närfält';%  On axis eller offaxis
VOLTAGE='';
MIKAVSTAND=100; %Mätavstånd till mikrofonen

disp('Trolldist')


%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
   
   
   
   
    DATA=REC(:,1);
    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_IMP.bmp',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

xlim([f_start f_end]);
ylim([20 110]);
STRING=sprintf('Frequency [Hz] | Trolldist v. %2.2f Beta',Trollver);
STRINGY=sprintf('SPL [dB] @ %d cm ref 20 uPa, %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.bmp',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')

%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 y 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('Utvolymen på förstärkaren är inte tillräckligt hög för önskad utnivå')
            disp(MaxA)
            return
        end
       
        %Spela in signalen
        disp('Recording signal')
        pause(0.1);
        wavplay(PLAY(:)*[1 1],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 y
       
        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);
       
        legend(STRI);
        drawnow
       
        %Spara plot till fil
        STRING=sprintf('%s_DIST.bmp',ELEMENT);
        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]);
    xlabel('Frequency [Hz]');
    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('DIST_%s.bmp',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


[Nefilim]

[phon]

Trolldist 7.3

860 rader, börjar bli svettigt.

Fixar många små saker. Bättre plottar. Beräknar själv de flesta hemliga konstaterna nu. Bestäm bara sveptid T f_start samt f_end.

Finns fortfarande många förbättringsmöjligheter

Ojoj, börjar bli mycke' nu

Jag har precis fått 7.0 att gå hyfsat på jobbet, tänkte man kanske skulle mäta på ACE-subben från AudioPro? Två 6,5 tummare och ett basrör (är det väl).

Såg en bild från Lillkubmätning med nån mik på, tänkte man skulle mäta ungefär lika. Hade du en mik eller två? Nära, 50cm, 1m, en vid varje element eller ...? Sitter en inbyggd stärkare där på 150w, den har nån sorts AGC-funktion också har jag för mig. Undrar hur det yttrar sig, intressant.

[Svante]

Två mikrofoner...

[lilltroll]

Två mikrofoner...

... fungerar prima när arean är densamma. Annars är det lite knepigare.