%   RandPhase.m


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


global Signal Fs Bits
[Signal,Fs,Bits]=wavread('Data:IADS:Sounds.wav:826.wav');

RealRandFac=0.2;
ImagRandFac=0.2;
                                                                    

[NPoints,NChan]=size(Signal);
Pow2=floor(log2(NPoints));
NPoints=2^Pow2;
Signal=Signal(1:NPoints,:);
FFTSignal=fft(Signal);
RealSignal=real(FFTSignal);
ImagSignal=imag(FFTSignal);
PhaseSignal=angle(FFTSignal);
PFFTSignal=abs(FFTSignal.*conj(FFTSignal))./NPoints;

RandVec1=randperm(NPoints);
RandVec2=randperm(NPoints);
NRealRandPoints=floor(NPoints.*RealRandFac);
if NRealRandPoints>0
	RealRandVec1=RandVec1(1:NRealRandPoints);
	RealRandVec2=RandVec2(1:NRealRandPoints);
	RealSignal(RealRandVec1,:)=RealSignal(RealRandVec2,:);
end
NImagRandPoints=floor(NPoints.*ImagRandFac);
if NImagRandPoints>0
	subplot(2,1,1)
	plot(ImagSignal(:,1))
	ImagRandVec1=RandVec1(1:NImagRandPoints);
	ImagRandVec2=RandVec2(1:NImagRandPoints);
	ImagSignal(ImagRandVec1,:)=ImagSignal(ImagRandVec2,:);
end
FFTRandSignal=RealSignal+i.*ImagSignal;
RandSignal=abs(ifft(FFTRandSignal));
PFFTRandSignal=abs(FFTRandSignal.*conj(FFTRandSignal))./NPoints;

subplot(3,1,1)
plot(PFFTSignal(1:NPoints./2,:))
subplot(3,1,2)
plot(PFFTRandSignal(1:NPoints./2,:))
subplot(3,1,3)
plot(PFFTSignal(1:NPoints./2,:)-PFFTRandSignal(1:NPoints./2,:))
pause(.001)
% PhaseSignal=(unwrap(angle(FFTSignal)).*180)./pi;

sound(RandSignal,Fs)
