现代信号处理.docx
《现代信号处理.docx》由会员分享，可在线阅读，更多相关《现代信号处理.docx（10页珍藏版）》请在冰豆网上搜索。
现代信号处理
4.信号的函数表达式为：
，其中，
为一随时间变化的随机过程，
为经过390—410Hz带通滤波器后的高斯白噪声，
为高斯白噪声，采样频率为1kHz，采样时间为2.048s。
（1）利用现代信号处理的知识进行信号谱估计；
（2）利用现代信号处理知识进行信号的频率提取；
（3）分别利用Winner滤波和Kalman滤波进行去噪；
（4）利用Wigner-Ville分布分析信号的时频特性。
（1）：
利用现代信号处理的知识进行信号谱估计：
经典谱估计中两种主要的方法为直接法和间接法，其中间接法则先根据N个样本数据
的样本自相关函数
（4.1）
其中
，且
。
计算样本自相关函数的Fourier变换，得到功率谱
（4.2）
周期图方法估计的功率谱为有偏估计，可通过加窗来减少其偏差。
定义为
（4.3）
式中
（4.4）
式中，
是窗函数
的Fourier变换。
功率谱估计程序为：
clear
clc
closeallhidden
sf=1000;nfft=2048;
t=0:
1/1000:
2.047;
A=normrnd（0,1,1,2048）;
N=wgn（1,2048,1）;
f1=390;f2=410;
wc1=2*f1/sf;
wc2=2*f2/sf;
%归一化频率
f0=[0wc1-0.05wc1wc2wc2+0.051];
B=[001100];%设置带通和带阻
weigh=[111];%设置带通和带阻权重
b=remez（50,f0,B,weigh）;%传函分子
D=filter（b,1,N）;
y=sin（2*pi*t*100）+1.5*sin（2*pi*t*300）+A.*sin（2*pi*t*200）+D+N;
a（1,:
）=y;a（2,:
）=y.*sin（y）;
x=a（1,:
）;
y=a（2,:
）-a（1,:
）;
f=0:
sf/nfft:
sf/2-sf/nfft;
w=boxcar（nfft）;%加矩形窗
z=psd（y,nfft,sf,w,nfft/2）;
nn=1:
nfft/2;
plot（f（nn）,abs（z（nn）））;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
gridon;
图4.1功率谱估计结果图
（2）.信号频率的提取用离散傅立叶算法
离散傅立叶算法程序
clear
clc
closeallhidden
sf=1000;nfft=2048;
t=0:
1/1000:
2.047;
A=normrnd（0,1,1,2048）;
N=wgn（1,2048,1）;
f1=390;f2=410;
wc1=2*f1/sf;
wc2=2*f2/sf;
%归一化频率
f0=[0wc1-0.05wc1wc2wc2+0.051];
B=[001100];%设置带通和带阻
weigh=[111];%设置带通和带阻权重
b=remez（50,f0,B,weigh）;%传函分子
D=filter（b,1,N）;
y=sin（2*pi*t*100）+1.5*sin（2*pi*t*300）+A.*sin（2*pi*t*200）+D+N;
t2=（0:
nfft-1）/sf;
f=（0:
nfft-1）*sf/nfft;
y1=abs（fft（y））;
f=f（1:
nfft/2）;
y1=y1（1:
nfft/2）;
plot（t,y）;
title（'原始信号'）;
axis（[02.047-68]）;
plot（f,y1）;
title（'fft频率提取'）;
axis（[050001000]）;
xlabel（'f/Hz'）;
gridon;
图4.2原始信号时域图
图4.3信号频谱
（3）分别利用Winner滤波和Kalman滤波进行去噪；
clearall
closeall
M=100;%维纳滤波器阶数
sf=1000;nfft=2048;
L=nfft;
t=0:
1/1000:
2.047;
A=normrnd（0,1,1,2048）;
N=wgn（1,2048,1）;
f1=390;f2=410;
wc1=2*f1/sf;
wc2=2*f2/sf;
%归一化频率f0
f0=[0wc1-0.05wc1wc2wc2+0.051];
B=[001100];%设置带通和带阻
weigh=[111];%设置带通和带阻权重
b=remez（50,f0,B,weigh）;%传函分子
D=filter（b,1,N）;
y=sin（2*pi*t*100）+1.5*sin（2*pi*t*300）+A.*sin（2*pi*t*200）+D+N;
phixx=xcorr（y,y）;
fori=1:
M
forj=1:
M
Rxx（i,j）=phixx（i-j+L）;
end
end
s=sin（2*pi*t*100）+1.5*sin（2*pi*t*300）+A.*sin（2*pi*t*200）;
phixs=xcorr（y,s）;
fori=1:
M
rxs（i）=phixs（i+L）;
end
h1=（inv（Rxx））*rxs';
%获得理想FIR滤波器系数h1
AA=sin（2*pi*t*100）+1.5*sin（2*pi*t*300）+A.*sin（2*pi*t*200）;
fori=1:
M
h（i）=AA（i）;
end
%绘图比较估计滤波器与实际滤波器
figure
k=1:
M;
plot（k,h（k）,'r',k,h1（k）,'b'）;
title（'Idealh（n）&Calculatedh（n）'）;
legend（'Idealh（n）','Calculatedh（n）'）;
xlabel（'n'）;ylabel（'h（n）'）;
%比较理想输出与实际输出
v=D+N;
S=conv（h,v）;
SI
（1）=S
（1）;
LL1=sin（2*pi*t*100）+1.5*sin（2*pi*t*300）+A.*sin（2*pi*t*200）;
fori=2:
L
SI（i）=LL1（i）;
end
figure
k=1:
L;
plot（k,s（k）,'r',k,SI（k）,'b'）;
title（'s（n）VS.SI（n）'）;
legend（'s（n）','SI（n）',0）;
xlabel（'n'）;ylabel（'IdealOutput'）;
holdon
SR=conv（h1,y）;
figure
k=1:
L;
plot（k,s（k）,'r',k,SR（k）,'b'）;
title（'s（n）VS.SR（n）'）;
legend（'s（n）去噪前','SR（n）去噪后',0）;
xlabel（'n'）;ylabel（'ActualOutput'）;
图4.4Winner滤波去噪图
Kalman滤波程序
clear;
clc;
Fs=1000;
nfft=2048;
t1=0:
1/Fs:
2.047;
A=normrnd（0,1,1,2048）;
N=wgn（1,2048,2）;
f1=390;f2=410;
wc1=2*f1/Fs;
wc2=2*f2/Fs;
wc2=2*f2/sf;
%归一化频率f0
f0=[0wc1-0.05wc1wc2wc2+0.051];
B=[001100];%设置带通和带阻
weigh=[111];%设置带通和带阻权重
b=remez（50,f0,B,weigh）;%传函分子
D=filter（b,1,N）;
x=sin（2*pi*t1*100）+1.5*sin（2*pi*t1*300）+A.*sin（2*pi*t1*200）+D+N;
x1=sin（2*pi*t1*100）+1.5*sin（2*pi*t1*300）+A.*sin（2*pi*t1*200）;
a1=-1.352;a2=1.338;a3=-0.662;a4=0.240;
A=[-a1-a2-a3-a4;1000;0100;0010];%状态转移矩阵
H=[1000];%观测矩阵
Q=[1000;0000;0000;0000];%状态噪声方差
R=1;%观测噪声方差阵
X（:
1）=[x（4）;x（3）;x
（2）;x
（1）];
p（:
:
1）=[10000;0100;0010;0001];%一步预测误差方针
%开始滤波
fork=2:
nfft
p1（:
:
k）=A*p（:
:
k-1）*A'+Q;%p1（:
:
k）即是一步预测误差的自相关矩阵，它是4*4的矩阵，取不同的k值就构成了一个三维矩阵
K（:
k）=p1（:
:
k）*H'/（H*p1（:
:
k）*H'+R）;%K（:
:
k）是增益矩阵,对于固定的k值它是4*1矩阵，取不同的k值就是三维矩阵
X（:
k）=A*X（:
k-1）+K（:
k）*[x（k）-H*A*X（:
k-1）];%X（:
k）是估计值,4*1矩阵
p（:
:
k）=p1（:
:
k）-K（:
k）*H*p1（:
:
k）;%p（:
:
k）是估计误差的自相关矩阵，4*4矩阵的三维矩阵
end%结束一次滤波
%绘图
t=1:
nfft;
figure
（2）;
plot（t,x1,'k-',t,x,'r-',t,X（1,:
）,'b-.'）;
title（'卡曼滤波去噪'）
legend（'真实轨迹','观测样本','估计轨迹'）;
gridon;
图5Kalman滤波去噪图
（4）利用Wigner-Ville分布分析信号的时频特性
MATLAB程序
clear;
clc;
Fs=1000;
nfft=2049;
t1=0:
1/Fs:
2.048;
A=normrnd（0,1,1,2049）;
N=wgn（1,2049,2）;
f1=390;f2=410;
wc1=2*f1/Fs;
wc2=2*f2/Fs;
%归一化频率f0
f0=[0wc1-0.05wc1wc2wc2+0.051];
B=[001100];%设置带通或带阻，1为带通，0为带阻
weigh=[111];%设置通带和阻带的权重
b=remez（50,f0,B,weigh）;%传函分子
D=filter（b,1,N）;
x=sin（2*pi*t1*100）+1.5*sin（2*pi*t1*300）+A.*sin（2*pi*t1*200）+D+N;
figure（8）
tfrwv（x'）;
xlabel（'时间t'）;
ylabel（'频率f'）;
图6幅频特性图
升级会员 