伪随机函数参数正弦波去噪处理

未知信号参数估计的应用研究

摘要：对于已接收到的某信号进行分析，此信号中也包含了噪声信息，利用测出的信号去估计原信号的的参数，参数暂定为幅度A、角频率ω和相位φ。利用已有的算法和Matlab的函数，进行去噪处理，同时对比两种算法的去噪效果，同时对比原信号函数和处理后的的信号，从图形和数值角度对所估计参数进行评估，此外对所得数据进行统计，得出两种算法得出参数估计的误差

关键词：参数估计，去噪算法，平滑函数，卷积函数

研究背景

现代高科技战争是海、陆、空、天与电、光、声信息一体化的战争，雷达在其中扮演着揭开战争迷雾天眼的角色。雷达回波信号处理作为雷达系统中的重要一环，对于已接受到的回波信号进行分析处理，目的是削弱信号中的多于内容：滤出混杂的噪声和干扰；或者将信号变为容易处理、传输、分析和识别的形式。雷达信号处理是指对观测到的信号进行分析、变换、综合等处理，以达到抑制干扰、杂波等非期望信号，增强有用信号，并估计有用信号的特征参数，或是将信号变成某种更符合要求的形式。随着微电子技术的迅速发展，信号处理的方式也从早期的模拟域发展到几乎都采用数字域。数字信号处理以数字或符号序列表示信号，用数值计算的方法完成对信号的各种处理。

处理原理

假设雷达发射信号为参数随机的正弦信号S(t)，其幅度A、角频率ω和相位φ未知，接收到雷达回波信号围为某未知信号Y(t)，同时近似将其看作为原信号S(t)和噪声N(t)的线性叠加。对回波信号进行有限个采样序列估计，估算出原信号的幅度、角频率和相位参数。下文将采用滑动平均法、卷积法和FIR滤波器处理噪声。

1.滑动平均法（moving average）

参考《数字信号处理》-胡广书中滑动平均法（moving average）也叫做移动平均法、平均法、移动平均值滤波法等等，是一种时间域思想上的信号光滑方法。算法思路为，将该点附近的采样点做算数平均，作为这个点光滑后的值。一般窗口为对称窗口，防止出现相位偏差。窗口一般为奇数。以3点平均（窗口长度为3）公式为例，原数据为x，平滑后的数据为y：y(n)=1/3∗(x(n−1)+x(n)+x(n+1))

2.卷积公式平滑法

根据以上的平滑公式，就相当于一个对x和向量[1/3 1/3 1/3]做卷积。对以上两个平滑法进行对比，结果如下图所示：

中间部分两者完全一致，但是两端有所差别。主要是因为，Movmean ()函数在处理边缘时，采用减小窗口的方式，而Conv ()相当于在两端补零。所以如何处理边缘也是值得注意的。

3.FIR滤波法实现光滑去噪

FIR滤波器也叫做 有限冲击响应滤波器。和它相对的是IIR 无限冲击响应滤波器。在Matlab里，有一张图可以比较直观。

对于FIR滤波器，用到的数据只有输入项b，输出项为1。但是对于IIR滤波器，不仅有输入项b（分子），还有输出项a（分母）。

对于FIR滤波器，由于只有分子b，所以可以用卷积Conv()函数去进行滤波。但是对于IIR，就不能用卷积函数去计算。一般认为，IIR函数在相同阶数下，滤波效果要比FIR函数要好，但是有可能出现发散问题。

比如前面的移动平均滤波中的1.3章所介绍的，b=[1/3,1/3,1/3]，a=1。就属于一个典型的FIR滤波器。其中，Conv(x,b)卷积方法相当于无相位偏移的中心滤波，

y(n)=1/3∗(x(n−1)+x(n)+x(n+1))

而filter(b,a,x)函数相当于有相位偏移的向后滤波

y(n)=1/3∗(x(n)+x(n+1)+x(n+2))

事实上，通过下面的等式，我们还可以构建一个等价的IIR滤波器：

正因为两个滤波器等价，其对信号的平滑处理效果也是相等的，效果图如下：

但是对比来看，FIR滤波对该信号的处理会比等价的IIR滤波更加接近原函数。

4.利用调制解调去噪。

调制：调制就是使一个信号（如光、高频电磁振荡等）的某些参数（如振幅、频率等）按照另一个欲传输的信号（如声音、图像等）的特点变化的过程。

解调：是调制的逆过程，作用是从已调信号中取出原来的调制信号。（解调则是将基带信号从载波中提取出来以便预定的接收者（信宿）处理和理解的过程。）

基带信号：要传递的信息，又称调制信号。

已调信号：调制波（信息）和载波叠加在一起就形成了已调信号。

载波：未受调制的周期性振荡信号称为载波，载波可以是正弦波，也可以是非正弦波（如周期性脉冲序列）。

一般要求正弦载波的频率远远高于调制信号的带宽，也就是常说的奈奎斯特采样频率，否则会发生混叠，使传输信号失真。

所得结果

对于未知参数的正弦函数信号，对其幅度，角频率和相位进行估计，在Matlab中模拟生成雷达回波信号，由正弦函数与高斯白噪声叠加。利用若干去噪方法对回波信号进行处理，对其分析得到参数估计得到其幅度，角频率和相位。

Step1

利用randi函数生成原信号参数，确保每次生成的源信号参数“随机”，但是可追踪变量求得此参数，以便后续的数值比较，再利用Matlab中的awgn( )函数进行高斯白噪声加扰，得到模拟雷达回波的信号Y(t)。如下图一所示:

图一.源信号（左）加扰信号（右）

上图中称为预处理，利用随机函数对参数做伪随机处理，也方便后续对比计算参数估计的准确性和误差大小。

其特点为每次运行后所得正弦函数参数不同，所得源信号不同，同时兼容通过追踪参数变量的值确定估算的准确标准

Step2

利用滑动平均法或叫移动评价值滤波法中的卷积函数Conv和Matlab自带Movmean函数对回波信号进行去噪处理,对比结果如下图二所示

图二，conv函数（左）和mov函数（右）

也正如前面原理部分所讲，movmean()函数在处理边缘时，采用减小窗口的方式，而conv()相当于在两端补零，在边缘处理时两个函数有所区别。

Step3

利用FIR滤波法和等价的IIR滤波法对其进行去噪处理，等价于前面所讲的a,b参数已知的中心滤波。对接收的雷达回波处理对比如图三所示：

图三.FIR和IIR滤波效果对比

对比以上的信号处理图，一般认为，IIR函数在相同阶数下，滤波效果要比FIR函数要好，但是有可能出现发散问题。前面的移动平均滤波中所介绍的，b=[1/3,1/3,1/3]，a=1。就属于一个典型的FIR滤波器。其中，conv(x,b)卷积方法相当于无相位偏移的中心滤波，而filter(b,a,x)函数相当于有相位偏移的向后滤波。

Step4

调剂与解调信号的调制简单来说，就是用基带信号去控制（调制）单频载波信号的某个或某几个参数，从而使得调制后的信号中嵌入了希望传递的信息。这个载波往往是一个高频信号，便于传输。解调则是调制的逆过程。但是更重要的可以在图中直观看出信号的频率，在解调后，频谱上高斯白噪声也被滤去。

图中可清晰看到频谱的搬移和复原，也能更准估算信号的频率约为7hz。

Step5

对相位处理，利用解调后的频谱相位图进行对比分析。

数据处理及结论

通过追踪Matlab中工作区中参数变化，同时再利用Findpeak( )函数对图像进行数据获取，对平滑处理后所的图像进行数据化表示，将所得数据导出求出其期望值和方差，与生成源函数进行比较，对比出去噪效果比较好且更趋向于原函数的算法。

原计划利用处理后的图像进行计算，但因其截断误差较大，后面又采取了老师所提出的调制解调方法，对信号进行处理，所得滤波的效果十分出色。但因采取的是SSB滤波，在幅值上所得信号仅为源信号的1/2，故将所得信号系数前乘常数2.在对其进行分析处理。进行三组对比

幅度	类别
3	源信号
2.945594154	解调
3.707694974	IIR
2.881516916	FIR
2.690548912	卷积
2.689523566	滑动

增加了观测时间取更多的值计算

对比可知相干解调对幅度估算最为实用，相对其他算法误差较小，同时还能得出频率和相位信息，是处理此类信号较为合适的的一种处理方法，对雷达回波信号进行傅里叶变化，在频域中获取信号所携带信息和参数特征

附录：Matlab代码

%%  清内存 关闭窗口 准备工作
clear 
close all
clc
 
%%  signal
% 三要素（ij均为整数）
A=randi([i  j]);              %amplify  幅度
f1=randi([i  j]);              %Hz  频率
w=2*pi*f1;                     %rad/s  角频率
m=randi([i j]);
p=pi/m;                       %rad  相位

%采样
T=time(整数);                %s        %观测时间
fs=20*f1;            %Hz       %采样频率
d=1/fs;             %s        %采样间隔
 
%正弦信号
t=0:d:T;                 %离散时间t
s1=A*sin(w*t+p);              %正弦信号
figure(1)
plot(t,s1,'r-')
 

%加噪白噪声
q1 = awgn(s1,10);
figure(2)
plot(t,q1)
 

%去噪后对比（利用移动平均法）
N_window = 5;%窗口长度(最好为奇数）
B1 = movmean(q1,N_window);%利用matlab自带函数
[pks1,locs1] = findpeaks(B1);
F_average = 1/N_window*ones(1,N_window);%构建卷积核
B2 = conv(q1,F_average,'same');%利用卷积的方法计算
[pks2,locs2] = findpeaks(B2);
figure(3)
plot(t,B1,'b--',t,B2,'r:')
legend('movmean','conv')
% FIR和IIR滤波
figure(4)
b=[1/3,1/3,1/3];a=1;
F1 = filter(b,a,q1);
[pks3,locs3] = findpeaks(F1);
b=1/3*[1,0,0,0,-1];a=[1,-1];
F2 = filter(b,a,q1);
[pks4,locs4] = findpeaks(F2);
plot(t,s1,'-r',t,F1,'-g',t,F2)
legend('原函数','FIR滤波','等价的IIR滤波')

%提取还原后和去噪后
figure(5)
plot(t,B1,'-g',t,s1,'-r');

legend('MOV去噪后','原函数')
figure(6)
plot(t,B2,'-b',t,s1,'-r');
set(gca,'ytick',-A:0.1:A);
grid on;    % 开启网格
legend('COV去噪后','原函数')
figure(7)
plot(t,s1,t,B1,t,F1,'-g')
legend('原函数','MOV光滑','FIR滤波')
% 调制解调
fm=f1;fc=40;
am=A;
Fs=fs;     %采样频率Fs，载波频率fc，信号频率fm
wc=2*pi*fc;
wm=2*pi*fm;
N=300;
n=0:N-1;
t=n/Fs;             %时间序列
f=n*Fs/N;
 
%基带信号时域
s2=am*sin(wm*t+p);
sm=awgn(s2,10);
S=fft(s2,300);%300点的fft
phase2 = angle(S);            % 相位谱（rad）
figure(8);
subplot(311);
plot(t,sm);
title('基带信号');
xlabel('t');
%axis([0 1 -2 2]);

 
%基带信号频域
S=fft(sm,300);%300点的fft
phase = angle(S);            % 相位谱（rad）
SG=abs(S);
subplot(312);
plot(f(1:N/2),SG(1:N/2));      %SSB信号频域波形
xlabel('Frequency(HZ)');
title('基带信号频域波形 ');
 

subplot(313);
plot(phase2, 'linewidth', 1);             %解调后的频域波形
title('源信号的频谱图');
xlabel('Frequency(HZ)');
grid on;
%SSB调制信号时域
s=modulate(sm,fc,Fs,'amssb');       %对调制信号进行调制
S=fft(s,300);
SG=abs(S);
figure(9);
subplot(211);
plot(t,s);                          %SSB信号时域波形
title('SSB信号时域波形 ');        
xlabel('t');
 
subplot(212);
plot(f(1:N/2),SG(1:N/2));            %SSB信号频域波形
xlabel('Frequency(HZ)');
title('SSB信号频域波形 ');
grid on;
%解调
fc=40;%载波频率fc
N=300;
n=0:N-1;
t=n/Fs;      %时间序列
f=n*Fs/N;
sm=am*cos(wm*t);
s=modulate(sm,fc,Fs,'amssb');  
sd=demod(s,fc,Fs,'amssb');         %对SSB信号进行解调
[pks5,locs5] = findpeaks(sd);
SD=fft(sd,300);
phase1 = angle(SD);                % 相位谱（rad）
SDG=abs(SD);
figure(10);
subplot(3,1,1);
plot(t,2*sd);                              %解调后的时域波形(特此乘2使得幅度相同)
title('解调后的时域波形');
xlabel('t');
grid on;
subplot(3,1,2);
plot(f(1:N/2),SDG(1:N/2));             %解调后的频域波形
title('解调后的频域波形');
xlabel('Frequency(HZ)');
grid on;
subplot(3,1,3);
plot(phase1, 'linewidth', 1);             %解调后的频域波形
title('解调得到的频谱图');
xlabel('Frequency(HZ)');
grid on;