gurencai 发表于 2009-3-27 13:06

产生标准高斯白噪声的算法

目前商业软件中产生的高斯白噪声序列都不是严格意义上的高斯白噪声，其均值，方差及正太分布性都不标准，此外白色性能也很差。本人现在在研究产生标准高斯白噪声的算法，文献中的大部分问题都已解决，但效果不是很好，希望有人能够帮忙或一起讨论。谢谢
QQ：504815336

gurencai 发表于 2009-3-31 14:31

白色化算法

http://fmn030.xnimg.cn/fmn030/pic001/20090331/14/30/large_4MJP_1790l241025.jpg

gurencai 发表于 2009-3-31 14:31

大家看看我的程序有什么问题没？

program standard_gauss_noise   !双随机交换法白色化高斯随机序列
   implicit none
   integer,parameter::N=511,maxN=2000
   real,parameter::kesi=0.1
   integer::i,j,r1,r2,k,a,b
   integer::p(2)
   real(kind=8)::seed_r,sum1,sum2
   real(kind=8)::x(0:N,0:maxN),r(0:N-1,0:maxN),s(0:maxN)
   a=0
   b=N
   seed_r=5.0
   do i=0,N
   !读取数据
      open(unit=10,file='gauss.txt')
      read(10,"(f12.8)") x(i,0)
   enddo

   do i=0,maxN-1
      sum2=0 !初始化
      s(i)=0.0
   if(i==0)   then !计算s(0)
      do k=0,N-1
      sum1=0.0
            do j=0,N-k-1
               sum1=x(j,i)*x(j+k,i)+sum1
            enddo
            r(k,i)=sum1/N
   !open(unit=122,file='r.txt')
   !write(122,*) r(k,i)
         enddo
         
         do k=1,N-1
            sum2=r(k,i)**2+sum2
         enddo
         s(i)=sum2
   write(*,*)s(i),i
   else
         !第一步：双随机交换产生新的随机序列
         call   nrabs(seed_r,a,b,p)
      r1=p(1)
      r2=p(2)
   !write(*,*) r1,r2
         do j=0,N
            x(j,i)=x(j,i-1)
         enddo
         call exchange(x(r1,i),x(r2,i))
      !第二步：计算序列的自相关函数
         do k=0,N-1
      sum1=0
            do j=0,N-k-1
               sum1=x(j,i)*x(j+k,i)+sum1
            enddo
            r(k,i)=sum1/N
   
         enddo
      !第三步：计算平方和
         do k=1,N-1
            sum2=r(k,i)**2+sum2
         enddo
         s(i)=sum2

   !第四步：若s(i)<kesi,则停止运算
         if(s(i)<kesi) then
         do j=0,N
         open(unit=225,file='D:\Download\matlab65\work\standerd_noise.txt')
         write(225,"(f12.8)") x(j,i)
      enddo
      write(*,*)i
         exit    !停止运算
         endif
         !第五步：若s(i)>=s(i-1),则双随机交换失败，放弃此次交换
         if(s(i)<s(i-1)) then
   open(unit=111,file='D:\Download\matlab65\work\s.txt')
         write(111,"(f12.8)")s(i)
         else
         call exchange(x(r1,i),x(r2,i)) !把原来换过的数据交换回来
   !write(*,*) i
         endif
   endif
   enddo
   stop
   end   program standard_gauss_noise
   

    subroutineexchange(a,b) !交换
    implicit none
    real(kind=8)::a,b,t
    t=a
    a=b
    b=t
    return
    end   


    subroutine nrabs(r,a,b,p) !产生两个随机整数
    integer,parameter::n=2
    integer p(n)
    integer a,b
real(kind=8)::r
    s=b-a+1.0
    k=log(s-0.5)/log(2.0)+1
    l=1
    do 10 i=1,k
10l=2*l
    k=1
    s=4.0*l
    i=1
20if((i<=l).and.(k<=n)) then
    r=r+r+r+r+r
m=r/s
r=r-m*s
j=a+r/4.0
if(j<=b) then
p(k)=j
k=k+1
end if
i=i+1
goto 20
end if
return
end

gurencai 发表于 2009-3-31 14:33

怎么这么久也没人回复啊？

无水1324 发表于 2009-4-1 15:50

回复 地板 gurencai 的帖子

这个问题不懂，你去Matlab区看看，怎么样？

octopussheng 发表于 2009-4-2 15:00

搜索论坛，有用matlab做高斯白噪声的

gurencai 发表于 2009-4-9 20:22

回复 板凳 gurencai 的帖子

顶上去。:@( :@( :@( :@( :@(
页: [1]
查看完整版本: 产生标准高斯白噪声的算法