randomsignal.m

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

clear
numberOfPoints = 1024;                           % 采样点数
fs = 44100;                                      % 采样频率  
fsOfSignal = 500;                                % 原信号的估计频率
ch = 1;                                          % 声音采样声道数

maxNumber = (1/fs) * (numberOfPoints - 1);
maxNumber2 = (1/fs) * (numberOfPoints/2 - 1);
x  = 0 : 1/fs : maxNumber;
x2 = 0 : 1/fs : maxNumber2;                      % 产生时域横轴


% RANDOM SIGNAL

% wave = wavrecord(numberOfPoints * 2,fs,ch,'double');
% for i = 1:1024
%     waveSignal(i) = wave(i+1024);
% end                                            % 声音信号采样

waveSignal = sin(fsOfSignal * 2 * pi * x);       % 用来模拟的周期正弦信号

ranSignal  = awgn(waveSignal,-7);                % 为周期信号叠加高斯白噪声 信噪比 1:2


figure(1)                                        % 显示原周期信号和混合高斯白噪声后的信号

subplot(2,1,1);
plot(x,waveSignal,'b-');
grid on;
title('原周期信号')

subplot(2,1,2);
plot(x,ranSignal,'b-');
grid on;
title('混合信号')

haveSignal = 1 - textWgn(ranSignal,numberOfPoints/2,numberOfPoints)    % 检测混合信号中是否存在有用信号 1表示有 0表示没有

equRS   = equality(ranSignal,numberOfPoints)                           
mVarRS  = meanVar(ranSignal,numberOfPoints)                           
varRS   = variance(equRS,mVarRS)                                      % 计算混合信号的均值、均方值、方差

figure(2);                                                             % 显示混合信号的均值、方差、频谱、自相关函数、功率谱

subplot(4,1,1);
plot(x,equRS,'b-');
grid on;
title('混合信号均值')

subplot(4,1,2);
plot(x,varRS,'b-');
grid on;
title('混合信号方差')

frequencyRS = fft(ranSignal);                                         % 计算混合信号的1024点FFT变换

dfRS = (1/(maxNumber/(numberOfPoints-1)))/numberOfPoints;
nRS  = 0:numberOfPoints-1;
fRS  = nRS*dfRS;                                                      % 进行横坐标的频域变换

subplot(4,1,3);
plot(fRS,abs(frequencyRS(nRS+1))*2/numberOfPoints,'r-');              % 进行纵坐标的修正，画出频谱
grid on;
title('混合信号频谱');
axis([0,2.25e4,0,1])

selfRelationRS = relation(ranSignal,ranSignal,equRS,equRS,varRS,varRS,numberOfPoints/2,numberOfPoints);  % 计算自相关函数

powerRS = fft(selfRelationRS);                                        % 功率谱密度是自相关函数的傅立叶变换
pdfRS   = (1/(maxNumber/(numberOfPoints/2-1)))/(numberOfPoints/2);
pnRS    = 0:numberOfPoints/2-1;
pfRS    = pnRS*pdfRS;

subplot(4,1,4);
plot(pfRS,abs(powerRS(pnRS+1))*2/(numberOfPoints/2),'k-');            % 计算功率谱密度
grid on;
title('功率谱密度');
axis([0,1.25e4,0,0.2])

figure(3);                                                            % 显示经过滤波后的信号，与原信号进行对比

equWS   = equality(waveSignal,numberOfPoints)
mVarWS  = meanVar(waveSignal,numberOfPoints)
varWS   = variance(equWS,mVarWS)
selfRelationWS = relation(waveSignal,waveSignal,equWS,equWS,varWS,varWS,numberOfPoints/2,numberOfPoints);   % 计算原信号各项数字特征

[b,a] = butter(9,fsOfSignal / (fs/2),'low');                           % Butterworth 低通滤波器 3dB截止频率等于原信号估计频率
freqz(b,a,128,1000);                                                  % 显示滤波器的幅频特性和相频特性

signalRE = filter(b,a,ranSignal);                                     % 进行滤波，恢复原信号

equRE   = equality(signalRE,numberOfPoints)                        
mVarRE  = meanVar(signalRE,numberOfPoints)
varRE   = variance(equRE,mVarRE)
selfRelationRE = relation(signalRE,signalRE,equRE,equRE,varRE,varRE,numberOfPoints/2,numberOfPoints);     % 计算恢复信号的各项数字特性

figure(4)

subplot(2,1,1)
plot(x,waveSignal,'b-');
grid on;
title('原信号')

subplot(2,1,2);
plot(x,signalRE,'b-');
grid on;
title('经过滤波后的信号')

figure(5);                                                         % 显示原信号各项数字特征

subplot(5,1,1);
plot(x,equWS,'b-');
grid on;
title('原信号均值')

subplot(5,1,2);
plot(x,varWS,'b-');
grid on;
title('原信号方差')

frequencyWS = fft(waveSignal);

dfWS = (1/(maxNumber/(numberOfPoints-1)))/numberOfPoints;
nWS  = 0:numberOfPoints-1;
fWS  = nWS*dfWS;
selfRelationWS = relation(waveSignal,waveSignal,equWS,equWS,varWS,varWS,numberOfPoints/2,numberOfPoints);

subplot(5,1,3);
plot(x2,selfRelationWS)
grid on;
title('原信号自相关函数')

subplot(5,1,4);
plot(fWS,abs(frequencyWS(nWS+1))*2/numberOfPoints,'r-');
grid on;
title('原信号频谱');
axis([0,2.25e4,0,1]);

powerWS = fft(selfRelationWS);
pdfWS   = (1/(maxNumber/(numberOfPoints/2-1)))/(numberOfPoints/2);
pnWS    = 0:numberOfPoints/2-1;
pfWS    = pnWS*pdfWS;

subplot(5,1,5);
plot(pfWS,abs(powerWS(pnWS+1))*2/(numberOfPoints/2),'k-');
grid on;
title('原信号功率谱密度');
axis([0,1.25e4,0,1]);

figure(6);                                                       % 显示恢复信号的各项数字特征

subplot(5,1,1);
plot(x,equRE,'b-');
grid on;
title('经过滤波后信号的均值')

subplot(5,1,2);
plot(x,varRE,'b-');
grid on;
title('经过滤波后信号的方差')

frequencyRE = fft(signalRE);

dfRE = (1/(maxNumber/(numberOfPoints-1)))/numberOfPoints;
nRE  = 0:numberOfPoints-1;
fRE  = nRE*dfRE;
selfRelationRE = relation(signalRE,signalRE,equRE,equRE,varRE,varRE,numberOfPoints/2,numberOfPoints);

subplot(5,1,3);
plot(x2,selfRelationRE);
grid on;
title('经过滤波后信号的自相关函数')

subplot(5,1,4);
plot(fRE,abs(frequencyRE(nRE+1))*2/numberOfPoints,'r-');
grid on;
title('经过滤波后信号的频谱');
axis([0,2.25e4,0,0.5]);


powerRE = fft(selfRelationRE);
pdfRE   = (1/(maxNumber/(numberOfPoints/2-1)))/(numberOfPoints/2);
pnRE    = 0:numberOfPoints/2-1;
pfRE    = pnRE*pdfRE;

subplot(5,1,5);
plot(pfRE,abs(powerRE(pnRE+1))*2/(numberOfPoints/2),'k-');
grid on;
axis([0,1.25e4,0,2]);

title('经过滤波后信号的功率谱密度');
 
periodWS = 1/fsOfSignal                                            % 计算原信号周期
[maxFsRE,i] = max(abs(frequencyRE(nRE+1))*2/numberOfPoints);
periodRE = 1/((i/numberOfPoints)*fs)                               % 计算恢复信号周期