现代信号处理.docx
- 文档编号:6224739
- 上传时间:2023-01-04
- 格式:DOCX
- 页数:10
- 大小:361.93KB
现代信号处理.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幅频特性图
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 信号 处理