弯弓射大雕 发表于 2007-9-15 10:31

看看我的滤波程序有啥问题,该怎样改进

我是用fortran编的,但是滤波结果与原始数据有较大差别,请高手指点,谢谢了
C       NUMERICAL FILTER 滤波
      
      PARAMETER(MN=5000000)
      DIMENSION AD(MN),DD(MN),IO(MN),DD2(MN)
      CHARACTER FINAME*15,FNA*15,SQUASH*40
      REAL*4    REAL(MN),IMAG(MN)
      OPEN(2,FILE=SQUASH(FINAME//'.800'))   
C      OPEN(12,FILE='AD.DAT')
      
      READ(2,*)(AD(I),I=1,NDH*NTOTAL) !读入数据
*****补零使数据长度为2的N次幂**************      
      MI=1
222   MNT=2**MI
      IF(MNT.LE.NTOTAL) THEN
          MI=MI+1
          GO TO 222
      ELSE      
          DO J=1,NDH
             DO I=NTOTAL+1,MNT
                K1=(I-1)*NDH+J
                K2=(I-2)*NDH+J
                AD(K1)=AD(K2)
             END DO   
         END DO
         NTOTAL=MNT
      END IF
          
      FN=1/((SCAN)*2.0)
       
      NFREQ=NTOTAL/FN
      DF=2./NFREQ
*****开始滤波******
         IFLT=1 !低通
      IF(IFLT.EQ.1) THEN
         WRITE(*,'(A\)')' ENTER HIGH CUT OFF FREQ.(Hz) : '
C         READ(*,*) FHI   
            FHI=0.02
         NHI=FHI/DF
      END IF
      IF(IFLT.EQ.2) THEN !高通
         WRITE(*,'(A\)')' ENTER LOW CUT OFF FREQ.(Hz) : '
         READ(*,*) FLO
         NLO=FLO/DF
      END IF
      IF(IFLT.EQ.3) THEN !带通
         WRITE(*,'(A\)')' ENTER HIGH CUT OFF FREQ.(Hz) : '
         READ(*,*) FHI
         NHI=FHI/DF
         WRITE(*,'(A\)')' ENTER LOW CUT OFF FREQ.(Hz) : '
         READ(*,*) FLO
         FLO=FLO/DF
      END IF
                  
****************************************************************         
      
      DO J=1,NDH
         DO I=1,NTOTAL
            K=(I-1)*NDH+J
            REAL(I)=AD(K)
            IMAG(I)=0.0
         END DO
         CALL FFT842(0,NTOTAL,REAL,IMAG) !快速fft
         IF (IFLT.EQ.1) THEN
            DO I=NHI+1,NTOTAL
               REAL(I)=0.0
               IMAG(I)=0.0
            END DO
         END IF
         IF (IFLT.EQ.2) THEN
            DO I=1,NLO
               REAL(I)=0.0
               IMAG(I)=0.0
            END DO
         END IF
         IF (IFLT.EQ.3) THEN
            DO I=NHI+1,NTOTAL
               REAL(I)=0.0
               IMAG(I)=0.0
            END DO
            DO I=1,NLO
               REAL(I)=0.0
               IMAG(I)=0.0
            END DO
         END IF
         CALL FFT842(1,NTOTAL,REAL,IMAG) !fft反变换
         DO I=1,NTOTAL
            K=(I-1)*NDH+J
            DD(K)=REAL(I)
            
         END DO
      END DO

弯弓射大雕 发表于 2007-9-15 10:37

这是滤波结果与原始数据的比较

弯弓射大雕 发表于 2007-9-15 10:39

老师说,滤波后信号失真了,我看也是滤波后与原来信号差别较大,请高手指点怎样改进呀
页: [1]
查看完整版本: 看看我的滤波程序有啥问题,该怎样改进