colorednoise.m

用matlab 实现强噪声中弱信号的提取

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('功率谱密度');