186 lines
		
	
	
		
			5.2 KiB
		
	
	
	
		
			Fortran
		
	
	
	
	
	
		
		
			
		
	
	
			186 lines
		
	
	
		
			5.2 KiB
		
	
	
	
		
			Fortran
		
	
	
	
	
	
|   | program fsk4sim
 | ||
|  | 
 | ||
|  |   parameter (ND=60)                      !Data symbols: LDPC (120,60), r=1/2
 | ||
|  |   parameter (NN=ND)                      !Total symbols (60)
 | ||
|  |   parameter (NSPS=57600)                 !Samples per symbol at 12000 sps
 | ||
|  |   parameter (NZ=NSPS*NN)                 !Samples in waveform (3456000)
 | ||
|  | 
 | ||
|  |   character*8 arg
 | ||
|  |   complex c(0:NZ-1)                      !Complex waveform
 | ||
|  |   complex cr(0:NZ-1)
 | ||
|  |   complex cs(NSPS,NN)
 | ||
|  |   complex cps(0:3)
 | ||
|  |   complex ct(0:2*NN-1)
 | ||
|  |   complex z,w,zsum
 | ||
|  |   real r(0:NZ-1)
 | ||
|  |   real s(NSPS,NN)
 | ||
|  |   real savg(NSPS)
 | ||
|  |   real tmp(NN)                          !For generating random data
 | ||
|  |   real xnoise(0:NZ-1)                   !Generated random noise
 | ||
|  |   real ps(0:3)
 | ||
|  |   integer id(NN)                         !Encoded 2-bit data (values 0-3)
 | ||
|  |   integer id2(NN)                        !Recovered data
 | ||
|  |   equivalence (r,cr)
 | ||
|  | 
 | ||
|  |   nnn=0
 | ||
|  |   nargs=iargc()
 | ||
|  |   if(nargs.ne.6) then
 | ||
|  |      print*,'Usage: fsk8sim f0 delay(ms) fspread(Hz) nts iters snr(dB)'
 | ||
|  |      go to 999
 | ||
|  |   endif
 | ||
|  |   call getarg(1,arg)
 | ||
|  |   read(arg,*) f0                        !Low tone frequency
 | ||
|  |   call getarg(2,arg)
 | ||
|  |   read(arg,*) delay
 | ||
|  |   call getarg(3,arg)
 | ||
|  |   read(arg,*) fspread
 | ||
|  |   call getarg(4,arg)
 | ||
|  |   read(arg,*) nts
 | ||
|  |   call getarg(5,arg)
 | ||
|  |   read(arg,*) iters
 | ||
|  |   call getarg(6,arg)
 | ||
|  |   read(arg,*) snrdb
 | ||
|  | 
 | ||
|  |   twopi=8.d0*atan(1.d0)
 | ||
|  |   fs=12000.d0
 | ||
|  |   dt=1.0/fs
 | ||
|  |   ts=NSPS*dt
 | ||
|  |   baud=1.d0/ts
 | ||
|  |   txt=NZ*dt
 | ||
|  |   bandwidth_ratio=2500.0/6000.0
 | ||
|  |   write(*,1000) baud,5*baud,txt,delay,fspread,nts
 | ||
|  | 1000 format('Baud:',f6.3,'  BW:',f5.1,'  TxT:',f5.1,'  Delay:',f5.2,   &
 | ||
|  |           '  fSpread:',f5.2,'  nts:',i3/)
 | ||
|  | 
 | ||
|  |   write(*,1004)
 | ||
|  | 1004 format('  SNR    Sym  Bit   SER     BER    Sym  Bit   SER     BER'/59('-'))
 | ||
|  | 
 | ||
|  |   isna=-25
 | ||
|  |   isnb=-40
 | ||
|  |   if(snrdb.ne.0.0) then
 | ||
|  |      isna=nint(snrdb)
 | ||
|  |      isnb=isna
 | ||
|  |   endif
 | ||
|  |   do isnr=isna,isnb,-1
 | ||
|  |      snrdb=isnr
 | ||
|  |      sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
 | ||
|  |      if(snrdb.gt.90.0) sig=1.0
 | ||
|  |      nhard1=0
 | ||
|  |      nhard2=0
 | ||
|  |      nbit1=0
 | ||
|  |      nbit2=0
 | ||
|  |      nh2=0
 | ||
|  |      nb2=0
 | ||
|  |      do iter=1,iters
 | ||
|  |         nnn=nnn+1
 | ||
|  |         id=0
 | ||
|  |         call random_number(tmp)
 | ||
|  |         where(tmp.ge.0.25 .and. tmp.lt.0.50) id=1
 | ||
|  |         where(tmp.ge.0.50 .and. tmp.lt.0.75) id=2
 | ||
|  |         where(tmp.ge.0.75) id=3
 | ||
|  | 
 | ||
|  |         call genfsk4(id,f0,nts,c)                !Generate the 4-FSK waveform
 | ||
|  |         call watterson(c,delay,fspread)
 | ||
|  |         if(sig.ne.1.0) c=sig*c                   !Scale to requested SNR
 | ||
|  |         if(snrdb.lt.90) then
 | ||
|  |            do i=0,NZ-1                           !Generate gaussian noise
 | ||
|  |               xnoise(i)=gran()
 | ||
|  |            enddo
 | ||
|  |         endif
 | ||
|  |         r(0:NZ-1)=real(c(0:NZ-1)) + xnoise       !Add noise to signal
 | ||
|  |         
 | ||
|  |         call snr2_wsprlf(r,freq,snr2500,width,1)
 | ||
|  |         write(*,3001) freq,snr2500,width
 | ||
|  | 3001    format(40x,3f10.3)
 | ||
|  |         
 | ||
|  |         df=12000.0/(2*NSPS)
 | ||
|  | !        i0=nint(f0/df)
 | ||
|  | !        i0=nint((1500.0+freq)/df)
 | ||
|  |         i0=nint((f0+freq)/df)
 | ||
|  |         call spec4(r,cs,s,savg)
 | ||
|  | 
 | ||
|  |         do j=1,NN
 | ||
|  |            nlo=0
 | ||
|  |            nhi=0
 | ||
|  |            ps=s(i0:i0+6*nts:2*nts,j)
 | ||
|  |            cps=cs(i0:i0+6*nts:2*nts,j)
 | ||
|  |            if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1
 | ||
|  |            if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1
 | ||
|  |            id2(j)=2*nhi+nlo
 | ||
|  |            z=cps(id2(j))
 | ||
|  |            ct(j-1)=z
 | ||
|  |         enddo
 | ||
|  |         nh1=count(id.ne.id2)
 | ||
|  |         nb1=count(iand(id,1).ne.iand(id2,1)) + count(iand(id,2).ne.iand(id2,2))
 | ||
|  | 
 | ||
|  |         ct(NN:)=0.
 | ||
|  |         call four2a(ct,2*NN,1,-1,1)
 | ||
|  |         df2=baud/(2*NN)
 | ||
|  |         ct=cshift(ct,NN)
 | ||
|  |         ppmax=0.
 | ||
|  |         dfpk=0.
 | ||
|  |         do i=0,2*NN-1
 | ||
|  |            f=(i-NN)*df2
 | ||
|  |            pp=real(ct(i))**2 + aimag(ct(i))**2
 | ||
|  |            if(pp.gt.ppmax) then
 | ||
|  |               ppmax=pp
 | ||
|  |               dfpk=f
 | ||
|  |            endif
 | ||
|  |         enddo
 | ||
|  | 
 | ||
|  |         zsum=0.
 | ||
|  |         do j=1,NN
 | ||
|  |            phi=(j-1)*twopi*dfpk*ts
 | ||
|  |            w=cmplx(cos(phi),sin(phi))
 | ||
|  |            cps=cs(i0:i0+6*nts:2*nts,j)*conjg(w)
 | ||
|  |            z=cps(id2(j))
 | ||
|  |            ct(j)=z
 | ||
|  |            zsum=zsum+z
 | ||
|  |            write(12,1042) j,id(j),id2(j),20*ps,atan2(aimag(z),real(z)),  &
 | ||
|  |                 atan2(aimag(zsum),real(zsum)),zsum
 | ||
|  | 1042       format(3i2,6f8.3,2f8.1)
 | ||
|  |         enddo
 | ||
|  | 
 | ||
|  |         phi0=atan2(aimag(zsum),real(zsum))
 | ||
|  |         zsum=0.
 | ||
|  |         do j=1,NN
 | ||
|  |            phi=(j-1)*twopi*dfpk*ts + phi0
 | ||
|  |            w=cmplx(cos(phi),sin(phi))
 | ||
|  |            nlo=0
 | ||
|  |            nhi=0
 | ||
|  |            cps=cs(i0:i0+6*nts:2*nts,j)*conjg(w)
 | ||
|  |            ps=real(cps)
 | ||
|  |            if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1
 | ||
|  |            if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1
 | ||
|  |            id2(j)=2*nhi+nlo
 | ||
|  |            z=cps(id2(j))
 | ||
|  |            ct(j)=z
 | ||
|  |            zsum=zsum+z
 | ||
|  |         enddo
 | ||
|  |         
 | ||
|  |         nh2=count(id.ne.id2)
 | ||
|  |         nb2=count(iand(id,1).ne.iand(id2,1)) + count(iand(id,2).ne.iand(id2,2))
 | ||
|  |         nhard1=nhard1+nh1
 | ||
|  |         nhard2=nhard2+nh2
 | ||
|  |         nbit1=nbit1+nb1
 | ||
|  |         nbit2=nbit2+nb2
 | ||
|  | 
 | ||
|  |         fdiff=1500.0+freq - f0
 | ||
|  |         write(13,1040)  snrdb,snr2500,f0,fdiff,width,nh1,nb1,nh2,nb2
 | ||
|  | 1040    format(2f7.1,f9.2,f7.2,f6.1,2(i8,i6))
 | ||
|  | 40      continue
 | ||
|  |      enddo
 | ||
|  | 
 | ||
|  |      ser1=float(nhard1)/(NN*iters)
 | ||
|  |      ser2=float(nhard2)/(NN*iters)
 | ||
|  |      ber1=float(nbit1)/(2*NN*iters)
 | ||
|  |      ber2=float(nbit2)/(2*NN*iters)
 | ||
|  |      write(*,1050)  snrdb,nhard1,nbit1,ser1,ber1,nhard2,nbit2,ser2,ber2
 | ||
|  |      write(14,1050) snrdb,nhard1,nbit1,ser1,ber1,nhard2,nbit2,ser2,ber2
 | ||
|  | 1050 format(f6.1,2(2i5,2f8.4))
 | ||
|  |   enddo
 | ||
|  |   write(*,1060) NN*iters,2*NN*iters
 | ||
|  | 1060 format(59('-')/'Max:  ',2i5)
 | ||
|  | 
 | ||
|  | 999 end program fsk4sim
 |