Parzen窗matlab仿真实验一、Parzen窗基本原理Parzen窗估计法是一种具有坚实理论基础和优秀性能的非参数函数估计方法,它能够较好地描述多维数据的分布状态。其基本思想就是利用一定范围内各点密度的平均值对总体密度函数进行估计。一般而言,设x为d维空间中任意一点,N是所选择的样本总数,为了对x处的分布概率密度^xp进行估计,以x为中心作一个边长为h的超立方体V,则其体积为dVh,为计算落入V中的样本数k,构造一个函数使得11,,1,2,...,20,jjduu其他(1)并使u满足条件0u,且1udu,则落入体积V中的样本数为1NiNixxkh,则此处概率密度的估计值是:^111NiixxNxpVh(2)式(2)是Parzen窗估计法的基本公式,1uV称为窗函数,或核函数。在Parzen窗估计法的基本公式中,窗宽h是一个非常重要的参数。当样本数N有限时,h对估计的效果有着较大的影响。二、Matlab实现Parzen窗仿真实验为了编写Parzen窗,首先要确定一个窗函数。一般可以选择的窗函数有方窗、正态窗等。接下来将分别用方窗、正态窗实现Parzen窗法。2.1、方窗1,1,2,...(,)20jjidihxxjdkxxh若其他(3)其中,h为超立方体的棱长。因为编程中所用的信号是一维的正态信号,所以此处d的取值为1。2.1.1方窗Parzen仿真实验流程图生成均值为0,方差为1,长度为N的一维正态随机信号用u判断ix是否在以x为中心h为长度的对称区间内则落入x为中心的区间内的样本个数加一YN用(2)式求出^xp绘制估计概率密度函数曲线2.1.2方窗Parzen仿真实验mtalab代码r=normrnd(0,1,1,10000);f=r(1:1000);f=sort(f);lth=1;h=[0.25,1,4];V1=h(1);V2=h(2);V3=h(3);V=[V1,V2,V3];a=zeros(1000,3);p=zeros(1000,3);b=0;form=1:3forx=1:1000forn=1:lthifabs((r(:,n)-f(x))/h(:,m))=1/2q=1;elseq=0;endb=q+b;enda(x,m)=b;b=0;endendform=1:3forx=1:1000p(x,m)=1/(lth*V(m))*a(x,m);endsubplot(4,3,m)plot(f(1:1000),p(:,m))ylabel('N=1')endholdonlth=16;b=0;form=1:3forx=1:1000forn=1:lthifabs((r(:,n)-f(x))/h(:,m))=1/2q=1;elseq=0;endb=q+b;enda(x,m)=b;b=0;endendform=1:3forx=1:1000p(x,m)=1/(lth*V(m))*a(x,m);endsubplot(4,3,m+3)plot(f(1:1000),p(:,m))ylabel('N=16')endholdonlth=256;b=0;form=1:3forx=1:1000forn=1:lthifabs((r(:,n)-f(x))/h(:,m))=1/2q=1;elseq=0;endb=q+b;enda(x,m)=b;b=0;endendform=1:3forx=1:1000p(x,m)=1/(lth*V(m))*a(x,m);endsubplot(4,3,m+6)plot(f(1:1000),p(:,m))axis([-5,5,0,0.6]);ylabel('N=256')endlth=10000;b=0;form=1:3forx=1:1000forn=1:lthifabs((r(:,n)-f(x))/h(:,m))=1/2q=1;elseq=0;endb=q+b;enda(x,m)=b;b=0;endendform=1:3forx=1:1000p(x,m)=1/(lth*V(m))*a(x,m);endsubplot(4,3,m+9)plot(f(1:1000),p(:,m))axis([-5,5,0,0.5]);ylabel('N=10000')end2.1.3方窗Parzen仿真实验运行结果h=0.25h=1h=42.2正态窗正态窗函数为21111(,)exp22ikxxuVVu(4)概率密度的估计式为2^111111exp22NiNixxNxphh(5)式中1hNh。2.2.1正态Parzen窗仿真实验流程图2.2.2正态Parzen窗仿真实验mtalab代码r=normrnd(0,1,1,10000);f=r(1:1000);f=sort(f);lth=1;h=[0.25,1,4];a=zeros(1000,3);p=zeros(1000,3);b=0;form=1:3h1=h(m)/sqrt(lth);forx=1:1000forn=1:lthb=b+exp(((r(:,n)-f(x))/h1).^2/(-2))/sqrt(2*pi)/h1;enda(x,m)=b/lth;b=0;endendform=1:3forx=1:1000生成均值为0,方差为1,长度为N的一维正态随机信号根据h和求得h1并计算出通过核函数的值用(5)式求出^xp绘制估计概率密度函数曲线p(x,m)=a(x,m);endsubplot(4,3,m)plot(f(1:1000),p(:,m))ylabel('N=1')endholdonlth=16;h=[0.25,1,4];a=zeros(1000,3);p=zeros(1000,3);b=0;form=1:3h1=h(m)/sqrt(lth);forx=1:1000forn=1:lthb=b+exp(((r(:,n)-f(x))/h1).^2/-2)/sqrt(2*pi)/h1;enda(x,m)=b/lth;b=0;endendform=1:3forx=1:1000p(x,m)=a(x,m);endsubplot(4,3,m+3)plot(f(1:1000),p(:,m))ylabel('N=16')endholdonlth=256;h=[0.25,1,4];a=zeros(1000,3);p=zeros(1000,3);b=0;form=1:3h1=h(m)/sqrt(lth);forx=1:1000forn=1:lthb=b+exp(((r(:,n)-f(x))/h1).^2/-2)/sqrt(2*pi)/h1;enda(x,m)=b/lth;b=0;endendform=1:3forx=1:1000p(x,m)=a(x,m);endsubplot(4,3,m+6)plot(f(1:1000),p(:,m))axis([-5,5,0,0.6]);ylabel('N=256')endlth=10000;h=[0.25,1,4];a=zeros(1000,3);p=zeros(1000,3);b=0;form=1:3h1=h(m)/sqrt(lth);forx=1:1000forn=1:lthb=b+exp(((r(:,n)-f(x))/h1).^2/-2)/sqrt(2*pi)/h1;enda(x,m)=b/lth;b=0;endendform=1:3forx=1:1000p(x,m)=a(x,m);endsubplot(4,3,m+9)plot(f(1:1000),p(:,m))axis([-5,5,0,0.5]);ylabel('N=10000')end2.2.3正态Parzen窗仿真实验运行结果h=0.25h=1h=42.3结论从方窗和正态窗两个Parzen窗仿真实验可以看出,估计的概率密度函数与N和h的取值大小有密切的关系。当h=0.25,N=10000时,曲线起伏很大,噪声大,与真实的概率密度函数曲线相差较大。h=1,起伏减小,不过噪声依旧明显。h=4,曲线平坦。同时当N→∞时,估计曲线也逐渐逼近实际的概率密度曲线。对比方窗与正态窗的两个Parzen仿真实验也可以看出,正态窗函数的估计结果要好于方窗,且光滑度也较好。(byzzd)