% demonstrate power spectrum of simple channel noise
% calculated by integrating a sequence of alternating 1's and -1's
% at random times
%
npts=10000; %number of samples
pchange=0.03; %probability of a state change
z=double(rand(1,npts)<pchange); %0's and 1's
zeven=find(z==1);
z(zeven(2:2:end))=-1; %z is now alternately +1 and -1
x=cumsum(z); %openings and closings
%
figure;
set(gcf,'Position',[100 200 1100 600]);
subplot('position',[0.1 0.8 0.8 0.1])
plot(x);set(gca,'YLim',[-0.2 1.2]);
title('simple channel noise');
subplot('position',[0.1 0.1 0.8 0.6])
[S1,f]=mtspectrumc(x,setfield([],'tapers',[1 1]));
loglog(f,S1,'c-'); hold on;
[S5,f]=mtspectrumc(x,setfield([],'tapers',[3 5]));
loglog(f,S5,'m-'); hold on;
[S13,f]=mtspectrumc(x,setfield([],'tapers',[7 13]));
loglog(f,S13,'r-'); hold on;
Spred=pchange./(4*pchange^2+(2.*pi*f).^2);
loglog(f,Spred,'k-');
title('power spectrum by multitaper method');
legend('NW=1 k=1','NW=3 k=5','NW=7 k=13','analytic','Location','SouthWest');