function [audiomodsignal,audiospektrum,audiomodspektrum,audiorealspektrum,audiomodrealspektrum,...
	audioimagspektrum,audiomodimagspektrum,audiomag,audiomodmag,audiophasen,audiomodphasen,phasenrandvec]...
	=randphase(audiosignal,audiofs,audiobits,amount)
	
%   randphase

%   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.                                                                      



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

% FFT-Größe gerade oder ungerade?
rest=rem(fftsize,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);

% FFT berechnen
audiospektrum=fft(audiosignal);

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

% Magnitude der Originalfrequenzenberechnen
for l=2:fftsize
    z1=(audiorealspektrum(l))^2;
    z2=(audioimagspektrum(l))^2;
    audiomag(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





% Zufällige Phasen generieren 

%Sinusförmige Verzerrung
%sinusindexvec=(linspace(0,1,fftsize))';
%phasenrandomvec=2*pi*sin(5*sinusindexvec);

%konstante Verzerrung
%phasenrandomvec=ones(fftsize,1)*0.5;

% zufällige Phasen = Rauschen
%phasenrandomvec=2*pi*(rand(fftsize,1));

%linear abnehmend
%phasenrandomvec=flipud(meshgrid(0:2/fftsize:6*pi,1)');

%frequenzabhängig
f=(1:fftsize/2)';

phasenrandomvec=(10000*ones(fftsize/2,1))./f;
phasenrandomvec(1:7000)=0;
phasenrandomvec(10000:fftsize/2)=0;

%sinusförmig
%a=linspace(pi/2,3*pi/2,fftsize/2);
%phasenrandomvec=2*pi*sin(a);

%hyperbel
%a=linspace(1,10,fftsize/2);
%c=ones(1,fftsize/2)*0.5;
%b=(rand(1,fftsize/2)-c)*20*pi;
%gradient=1./a.^10;
%phasenrandomvec=b.*gradient;



phasenrandvec=zeros(fftsize,1);  
% für Frequenzen bis NQ-Frequenz mit Originalphasen gewichtet mitteln
if rest==0
   for o=2:fftsize/2
      phasenrandvec(o)=(audiophasen(o)+phasenrandomvec(o));
   end
else
   for o=2:round(fftsize/2-1)
      phasenrandvec(o)=(audiophasen(o)+phasenrandomvec(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))*audiomag(j); 
         audiomodrealspektrum(j,1)=sqrt(((audiomag(j))^2)-((audiomodimagspektrum(j,1))^2));
      elseif and(phasenrandvec(j)>(pi/2),phasenrandvec(j)<=pi)
         audiomodimagspektrum(j,1)=sin(pi-phasenrandvec(j))*audiomag(j); 
       	audiomodrealspektrum(j,1)=(-1)*(sqrt(((audiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      elseif and(phasenrandvec(j)>pi,phasenrandvec(j)<=(3*pi/2))
         audiomodimagspektrum(j,1)=(-1)*sin(phasenrandvec(j)-pi)*audiomag(j);
         audiomodrealspektrum(j,1)=(-1)*(sqrt(((audiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      else 
         audiomodimagspektrum(j,1)=(-1)*sin(2*pi-phasenrandvec(j))*audiomag(j);
         audiomodrealspektrum(j,1)=sqrt((audiomag(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))*audiomag(j); 
         audiomodrealspektrum(j,1)=sqrt(((audiomag(j))^2)-((audiomodimagspektrum(j,1))^2));
      elseif and(phasenrandvec(j)>(pi/2),phasenrandvec(j)<=pi)
         audiomodimagspektrum(j,1)=sin(pi-phasenrandvec(j))*audiomag(j); 
        	audiomodrealspektrum(j,1)=(-1)*(sqrt(((audiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      elseif and(phasenrandvec(j)>pi,phasenrandvec(j)<=(3*pi/2))
         audiomodimagspektrum(j,1)=(-1)*sin(phasenrandvec(j)-pi)*audiomag(j);
         audiomodrealspektrum(j,1)=(-1)*(sqrt(((audiomag(j))^2)-((audiomodimagspektrum(j,1))^2)));
      else 
         audiomodimagspektrum(j,1)=(-1)*sin(2*pi-phasenrandvec(j))*audiomag(j);
         audiomodrealspektrum(j,1)=sqrt((audiomag(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)=audiospektrum(1,1);

% Nyquistfrequenz kopieren
if rest==0
   audiomodspektrum(fftsize/2+1,1)=real(audiospektrum(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));


fftsize
return
