63 lines
		
	
	
		
			1.5 KiB
		
	
	
	
		
			Fortran
		
	
	
	
	
	
		
		
			
		
	
	
			63 lines
		
	
	
		
			1.5 KiB
		
	
	
	
		
			Fortran
		
	
	
	
	
	
| 
								 | 
							
								subroutine dopspread(c,fspread)
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								  parameter (NFFT=268800,NH=NFFT/2)
							 | 
						||
| 
								 | 
							
								  complex c(0:NFFT-1)
							 | 
						||
| 
								 | 
							
								  complex cspread(0:NFFT-1)
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								  df=12000.0/nfft
							 | 
						||
| 
								 | 
							
								  twopi=8*atan(1.0)
							 | 
						||
| 
								 | 
							
								  cspread(0)=1.0
							 | 
						||
| 
								 | 
							
								  cspread(NH)=0.
							 | 
						||
| 
								 | 
							
								  b=6.0                                     !Lorenzian 3/28 onward
							 | 
						||
| 
								 | 
							
								  do i=1,NH
							 | 
						||
| 
								 | 
							
								     f=i*df
							 | 
						||
| 
								 | 
							
								     x=b*f/fspread
							 | 
						||
| 
								 | 
							
								     z=0.
							 | 
						||
| 
								 | 
							
								     a=0.
							 | 
						||
| 
								 | 
							
								     if(x.lt.3.0) then                          !Cutoff beyond x=3
							 | 
						||
| 
								 | 
							
								        a=sqrt(1.111/(1.0+x*x)-0.1)             !Lorentzian
							 | 
						||
| 
								 | 
							
								        call random_number(r1)
							 | 
						||
| 
								 | 
							
								        phi1=twopi*r1
							 | 
						||
| 
								 | 
							
								        z=a*cmplx(cos(phi1),sin(phi1))
							 | 
						||
| 
								 | 
							
								     endif
							 | 
						||
| 
								 | 
							
								     cspread(i)=z
							 | 
						||
| 
								 | 
							
								     z=0.
							 | 
						||
| 
								 | 
							
								     if(x.lt.50.0) then
							 | 
						||
| 
								 | 
							
								        call random_number(r2)
							 | 
						||
| 
								 | 
							
								        phi2=twopi*r2
							 | 
						||
| 
								 | 
							
								        z=a*cmplx(cos(phi2),sin(phi2))
							 | 
						||
| 
								 | 
							
								     endif
							 | 
						||
| 
								 | 
							
								     cspread(NFFT-i)=z
							 | 
						||
| 
								 | 
							
								  enddo
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								  izh=fspread/df
							 | 
						||
| 
								 | 
							
								  do i=-izh,izh
							 | 
						||
| 
								 | 
							
								     f=i*df
							 | 
						||
| 
								 | 
							
								     j=i
							 | 
						||
| 
								 | 
							
								     if(j.lt.0) j=j+nfft
							 | 
						||
| 
								 | 
							
								     s=real(cspread(j))**2 + aimag(cspread(j))**2
							 | 
						||
| 
								 | 
							
								!          write(23,3000) f,s,cspread(j)
							 | 
						||
| 
								 | 
							
								!3000      format(f10.3,3f12.6)
							 | 
						||
| 
								 | 
							
								  enddo
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								  call four2a(cspread,NFFT,1,1,1)             !Transform to time domain
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								  sum=0.
							 | 
						||
| 
								 | 
							
								  do i=0,NFFT-1
							 | 
						||
| 
								 | 
							
								     p=real(cspread(i))**2 + aimag(cspread(i))**2
							 | 
						||
| 
								 | 
							
								     sum=sum+p
							 | 
						||
| 
								 | 
							
								  enddo
							 | 
						||
| 
								 | 
							
								  avep=sum/NFFT
							 | 
						||
| 
								 | 
							
								  fac=sqrt(1.0/avep)
							 | 
						||
| 
								 | 
							
								  cspread=fac*cspread                   !Normalize to constant avg power
							 | 
						||
| 
								 | 
							
								  c=cspread*c                           !Apply Rayleigh fading to c()
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								  do i=0,NFFT-1
							 | 
						||
| 
								 | 
							
								     p=real(cspread(i))**2 + aimag(cspread(i))**2
							 | 
						||
| 
								 | 
							
								!     write(24,3010) i,p,cspread(i)
							 | 
						||
| 
								 | 
							
								!3010 format(i8,3f12.6)
							 | 
						||
| 
								 | 
							
								  enddo
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								  return
							 | 
						||
| 
								 | 
							
								end subroutine dopspread
							 |