function [audiomag,kritaudiomag,audiomodmag]=ampmatch(filename,filepath,kritfilename,kritfilepath) 

%   ampmatch

%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      

cd (filepath);
[audio.signal,audio.fs,audio.bits]=wavread(filename);
audiosignal=audio.signal;
audiofs=audio.fs;
audiobits=audio.bits;

cd (kritfilepath);
[kritaudio.signal,kritaudio.fs,kritaudio.bits]=wavread(kritfilename);
kritaudiosignal=kritaudio.signal;
kritaudiofs=kritaudio.fs;
kritaudiobits=kritaudio.bits;


% FFT-Größe ermitteln
fftsizevec=size(audiosignal);
fftsize=fftsizevec(1);

kritfftsizevec=size(kritaudiosignal);
kritfftsize=kritfftsizevec(1);


% FFT-Größe gerade oder ungerade?
rest=rem(fftsize,2);
kritrest=rem(kritfftsize,2);

% Matritzen für modifiertes Spektrum und Signal und für modifizierten real- und imag-Teil und erstellen
audiomodrealspektrum=zeros(fftsize,1);
audiomodimagspektrum=zeros(fftsize,1);  
audiomodspektrum=zeros(fftsize,1);
audiomodsignal=zeros(fftsize,1);
audiophasen=zeros((fftsize),1);
audiomodphasen=zeros((fftsize),1);
audiomag=zeros(fftsize,1);
audiomodmag=zeros(fftsize,1);

kritaudiomodrealspektrum=zeros(fftsize,1);
kritaudiomodimagspektrum=zeros(fftsize,1);  
kritaudiomodspektrum=zeros(fftsize,1);
kritaudiomodsignal=zeros(fftsize,1);
kritaudiophasen=zeros((fftsize),1);
kritaudiomodphasen=zeros((fftsize),1);
kritaudiomag=zeros(fftsize,1);
kritaudiomodmag=zeros(fftsize,1);


% FFT berechnen
audiospektrum=fft(audiosignal);
kritaudiospektrum=fft(kritaudiosignal);

% Komplexe Zahlen in real und imag-Teil aufsplitten
audiorealspektrum=real(audiospektrum);
audioimagspektrum=imag(audiospektrum);
kritaudiorealspektrum=real(kritaudiospektrum);
kritaudioimagspektrum=imag(kritaudiospektrum);

% Magnitude der Originalfrequenzenberechnen
for l=2:fftsize
    z1=(audiorealspektrum(l))^2;
    z2=(audioimagspektrum(l))^2;
    audiomag(l,1)=sqrt(z1+z2);
end

for l=2:kritfftsize
    z1=(kritaudiorealspektrum(l))^2;
    z2=(kritaudioimagspektrum(l))^2;
    kritaudiomag(l,1)=sqrt(z1+z2);
end


% Originalphasen berechnen
for p=2:fftsize
      if and(audioimagspektrum(p)>=0,audiorealspektrum(p)>=0)
         if audiomag(p)>0
            audiophasen(p)=asin((audioimagspektrum(p))/audiomag(p));
         end
      elseif and(audioimagspektrum(p)<0,audiorealspektrum(p)>=0)
         if audiomag(p)>0
            audiophasen(p)=asin((audioimagspektrum(p))/audiomag(p))+2*pi;
         end
      elseif and(audioimagspektrum(p)>=0,audiorealspektrum(p)<0)
         if audiomag(p)>0
            audiophasen(p)=pi-asin((audioimagspektrum(p))/audiomag(p));
         end
      elseif and(audioimagspektrum(p)<0,audiorealspektrum(p)<0)
         if audiomag(p)>0
            audiophasen(p)=abs(asin((audioimagspektrum(p))/audiomag(p)))+pi;
         end
      end
end

for p=2:kritfftsize
      if and(kritaudioimagspektrum(p)>=0,kritaudiorealspektrum(p)>=0)
         if kritaudiomag(p)>0
            kritaudiophasen(p)=asin((kritaudioimagspektrum(p))/kritaudiomag(p));
         end
      elseif and(kritaudioimagspektrum(p)<0,kritaudiorealspektrum(p)>=0)
         if kritaudiomag(p)>0
            kritaudiophasen(p)=asin((kritaudioimagspektrum(p))/kritaudiomag(p))+2*pi;
         end
      elseif and(kritaudioimagspektrum(p)>=0,kritaudiorealspektrum(p)<0)
         if kritaudiomag(p)>0
            kritaudiophasen(p)=pi-asin((kritaudioimagspektrum(p))/kritaudiomag(p));
         end
      elseif and(kritaudioimagspektrum(p)<0,kritaudiorealspektrum(p)<0)
         if kritaudiomag(p)>0
            kritaudiophasen(p)=abs(asin((kritaudioimagspektrum(p))/kritaudiomag(p)))+pi;
         end
      end
end





% zufällige Phasen = Rauschen
%phasenrandomvec=2*pi*(rand(fftsize,1));
%phasenrandvec=zeros(fftsize,1);
   
   
% Phasen beibehalten
if rest==0
   for o=2:fftsize/2
      phasenrandvec(o)=audiophasen(o);
   end
else
   for o=2:round(fftsize/2-1)
      phasenrandvec(o)=audiophasen(o);
   end
end


% modifierten real- und iaudiomag-Teil bis NQ-Frequenz berechnen  
if rest==0
   for j=2:fftsize/2
      if phasenrandvec(j)<(pi/2)
         audiomodimagspektrum(j,1)=sin(phasenrandvec(j))*kritaudiomag(j); 
         audiomodrealspektrum(j,1)=sqrt(((kritaudiomag(j))^2)-((audiomodimagspektrum(j,1))^2));
      elseif and(phasenrandvec(j)>(pi/2),phasenrandvec(j)<=pi)
         audiomodimagspektrum(j,1)=sin(pi-phasenrandvec(j))*kritaudiomag(j); 
       	audiomodrealspektrum(j,1)=(-1)*(sqrt(((kritaudiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      elseif and(phasenrandvec(j)>pi,phasenrandvec(j)<=(3*pi/2))
         audiomodimagspektrum(j,1)=(-1)*sin(phasenrandvec(j)-pi)*kritaudiomag(j);
         audiomodrealspektrum(j,1)=(-1)*(sqrt(((kritaudiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      else 
         audiomodimagspektrum(j,1)=(-1)*sin(2*pi-phasenrandvec(j))*kritaudiomag(j);
         audiomodrealspektrum(j,1)=sqrt((kritaudiomag(j)^2)-((audiomodimagspektrum(j,1))^2));
     	end          
  	end   
else 
    for j=2:round(fftsize/2-1)
      if phasenrandvec(j)<(pi/2)
         audiomodimagspektrum(j,1)=sin(phasenrandvec(j))*kritaudiomag(j); 
         audiomodrealspektrum(j,1)=sqrt(((kritaudiomag(j))^2)-((audiomodimagspektrum(j,1))^2));
      elseif and(phasenrandvec(j)>(pi/2),phasenrandvec(j)<=pi)
         audiomodimagspektrum(j,1)=sin(pi-phasenrandvec(j))*kritaudiomag(j); 
        	audiomodrealspektrum(j,1)=(-1)*(sqrt(((kritaudiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      elseif and(phasenrandvec(j)>pi,phasenrandvec(j)<=(3*pi/2))
         audiomodimagspektrum(j,1)=(-1)*sin(phasenrandvec(j)-pi)*kritaudiomag(j);
         audiomodrealspektrum(j,1)=(-1)*(sqrt(((kritaudiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      else 
         audiomodimagspektrum(j,1)=(-1)*sin(2*pi-phasenrandvec(j))*kritaudiomag(j);
         audiomodrealspektrum(j,1)=sqrt((kritaudiomag(j)^2)-((audiomodimagspektrum(j,1))^2));
     	end          
  	end    
end

% modifierten real und imag-Teil bis NQ-Frequenz zusammenfügen
if rest==0
   for q=2:fftsize/2
      audiomodspektrum(q)=audiomodrealspektrum(q) + (audiomodimagspektrum(q))*sqrt(-1);
   end
else
   for q=2:round(fftsize/2)
      audiomodspektrum(q)=audiomodrealspektrum(q) + (audiomodimagspektrum(q))*sqrt(-1);
   end
end

% Baselevel kopieren
audiomodspektrum(1,1)=kritaudiospektrum(1,1);

% Nyquistfrequenz kopieren
if rest==0
   audiomodspektrum(fftsize/2+1,1)=real(kritaudiospektrum(fftsize/2+1,1));
end

% um conjugates ergänzen
if rest==0
   startcon=fftsize/2+2;
   for j=startcon:fftsize
   	audiomodspektrum(j,1)=conj(audiomodspektrum((fftsize-j+2),1));
	end
else
   startcon=round(fftsize/2);
   for j=startcon:fftsize
   	audiomodspektrum(j,1)=conj(audiomodspektrum((fftsize-j+2),1));
	end
end

% Modphasen berechnen
for p=2:fftsize
      if and(audiomodimagspektrum(p)>=0,audiomodrealspektrum(p)>=0)
         if audiomodmag(p)>0
            audiomodphasen(p)=asin((audiomodimagspektrum(p))/audiomodmag(p));
         end
      elseif and(audiomodimagspektrum(p)<0,audiomodrealspektrum(p)>=0)
         if audiomodmag(p)>0
            audiomodphasen(p)=asin((audiomodimagspektrum(p))/audiomodmag(p))+2*pi;
         end
      elseif and(audiomodimagspektrum(p)>=0,audiomodrealspektrum(p)<0)
         if audiomodmag(p)>0
            audiomodphasen(p)=pi-asin((audiomodimagspektrum(p))/audiomodmag(p));
         end
      elseif and(audiomodimagspektrum(p)<0,audiomodrealspektrum(p)<0)
         if audiomodmag(p)>0
            audiomodphasen(p)=asin((audiomodimagspektrum(p))/audiomodmag(p))+pi;
         end
      end
end

% Ergebnisse über in real und imag-Teil aufsplitten
audiomodrealspektrum=real(audiomodspektrum);
audiomodimagspektrum=imag(audiomodspektrum);

% Magnitude der modifierten Frequenzen berechnen
for l=1:fftsize
    zm1=(audiomodrealspektrum(l,1))^2;
    zm2=(audiomodimagspektrum(l,1))^2;
    audiomodmag(l,1)=sqrt(zm1+zm2);
end
 
% Inverse FFT berechnen, modifiziertes Signal erstellen 
audiomodsignal=real(ifft(audiomodspektrum));

%stereo-sound
audiomodsignal=[audiomodsignal audiomodsignal];
sound(audiomodsignal,audiofs);

% modifiziertes Signal abspeichern
cd (filepath);
targetfilename=strrep(filename, '.wav','match.wav');
wavwrite(audiomodsignal,audiofs,audiobits,targetfilename);
   
   
return
   
   
   