% script to make fixed amplitude-random phase noise with sinusoidal AM
clear all;
close all;

srate=44100;
dur=input('enter sound duration in seconds ');
f1=input('enter lower passband in Hz ');
f2=input('enter upper passband in Hz ');
phi=input('enter modulation rate in Hz ');
depth=input('enter modulation depth as % ');
wavname=input('enter name for WAV File to be generated! ','s');

t=0:1/srate:dur-1/srate;

% number of samples should be even
if mod(length(t),2)~=0,
        t=0:1/srate:dur;
end

npts=length(t);
fbin=srate/npts;

mag=ones(1,npts/2);
mag(1:floor(f1/fbin))=0;
mag(floor(f2/fbin):length(mag))=0;

phase=2*pi*rand(1,npts/2-1);
allphase=[0 phase 0 -1.*fliplr(phase)];
allmag=[mag fliplr(mag)];
rect=allmag.*exp(i.*allphase);
sig=real(ifft(rect));

% modulate the noise


mod = 1 + ((depth/100)*(sin(2*pi*phi*t)));

sig = mod .* sig;


% set rms value to 0.1
sig=sig/(10*std(sig));

figure(1), subplot(2,1,1),
plot(t,sig);
xlabel('Time (s)');
ylabel('Digital Amplitude');

spc=20*log10(abs(fft(sig)));
spc=spc(1:npts/2);
spc=spc-max(spc);
frqs=0:fbin:srate/2-fbin;
subplot(2,1,2), plot(frqs/1e3,spc);
set(gca,'xtick',[0:2:20]);
set(gca,'ylim',[-60 0]);
set(gca,'xlim',[0 srate/2e3]);
xlabel('Frequency (kHz)');
ylabel('Magnitude (dB)');

wavwrite(sig,srate,wavname);
sound(sig,srate)