numberOfPoints = 10000;                                        % 采樣點數
maxNumber = 1;                                                 % 采樣時間

numberOfPictures = 5;

x = 0 : ((maxNumber)/(numberOfPoints-1)) : maxNumber;
x2 = 0 : ((maxNumber)/(numberOfPoints/2-1)) : maxNumber;
i = 1:numberOfPoints;


% COLORED NOISE
sawSignal = sawtooth(x);       % 產生鋸齒波
coloredNoise = awgn(x,10,'measured');                          % 將鋸齒波疊加白噪聲產生色噪聲

figure(8)

plot(x,sawSignal)

figure(9);                                                     % 顯示色噪聲及各項數字特征

subplot(numberOfPictures,1,1);
plot(x,coloredNoise);
grid on;
% axis([0,50,-5,5]);
title('色噪聲');

isCN = textWgn(coloredNoise,numberOfPoints/2,numberOfPoints)  % 測試是否是白噪聲  1表示是 0表示不是

equCN  = equality(coloredNoise,numberOfPoints)
mVarCN = meanVar(coloredNoise,numberOfPoints)
varCN  = variance(equCN,mVarCN)                                % 計算白噪聲的均值、均方值、方差

subplot(numberOfPictures,1,2);
plot(x,varCN,'k-');
grid on;
title('方差');

frequencyCN = fft(coloredNoise);                              % 計算色噪聲的FFT變換

dfCN = (1/(maxNumber/(numberOfPoints-1)))/numberOfPoints;
nCN  = 0:numberOfPoints-1;
fCN  = nCN*dfCN;                                              % 將橫軸變換到頻域

subplot(numberOfPictures,1,3);
plot(fCN,abs(frequencyCN(nCN+1))*2/numberOfPoints,'r-');      % 修正縱軸，畫出白噪聲頻譜
grid on;
title('頻譜');
axis([0,10000,0,0.05]);

selfRelationCN = relation(coloredNoise,coloredNoise,equCN,equCN,varCN,varCN,numberOfPoints/2,numberOfPoints);

subplot(numberOfPictures,1,5);
plot(x2,selfRelationCN,'k-');
grid on;
%axis([0,50,-0.5,0.5]);
title('自相關函數');
 
powerCN = fft(selfRelationCN);                                      % 功率譜密度是自相關函數的傅立葉變換
pdfCN   = (1/(maxNumber/(numberOfPoints/2-1)))/(numberOfPoints/2);
pnCN    = 0:numberOfPoints/2-1;
pfCN    = pnCN*pdfCN;

subplot(numberOfPictures,1,4); 
plot(pfCN,abs(powerCN(pnCN+1))*2/(numberOfPoints/2),'k-');          % 畫出白噪聲的功率譜密度
grid on;
axis([0,numberOfPoints/2,-0.05,0.05]);
title('功率譜密度');