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

gurencai

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

gurencai

                               
登录/注册后可看大图

gurencai

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
   

    subroutine  exchange(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
10  l=2*l
    k=1
    s=4.0*l
    i=1
20  if((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

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

无水1324

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

octopussheng

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

gurencai

顶上去。:@( :@( :@( :@( :@(