Silverfrost Forums

Welcome to our forums

code to window averaging filter

4 Jul 2011 9:05 #8509

Hi there I have a two-column data files with 8192 rows. Can someone tell me how to average each 8 values as tried below? I managed to average first eight values but do not know how to call the succeeding 8 values. The whole purpose is to remove noise from pressure data by window averaging. Thanks!

program win_avg_filter
implicit none
integer::n,i
real,dimension(8192)::y
open (unit=3, FILE='ptdata.txt', stat=ierror)
ierror = 0
do i = 1,8192
read (3,*,stat=ierror) ydummy,y(i)  
if (ierror /=0) exit
enddo
close (unit=3)
do i = 1,8
total = total + total(i)
fvalue = total/8.

[/code]

7 Jul 2011 9:56 #8530

Quoted from Srabon Hi there I have a two-column data files with 8192 rows. Can someone tell me how to average each 8 values as tried below? I managed to average first eight values but do not know how to call the succeeding 8 values. The whole purpose is to remove noise from pressure data by window averaging. Thanks!

This is my version. I have assumed the window is always 8 values and cannot move past the ends of the data. The output array (z) therefore has 8192 - 7 = 8185 elements. I used the sum intrinsic to add up the element values.

Incidently, I have not tried to compile this yet, I just typed it in here. I hope it compiles :lol:

I provide a more optimsed version below.

program win_avg_filter
   implicit none

   integer, parameter :: WINSIZ = 8
   integer:: i,  ierror
   real ::y(8192),  z(8192-WINSIZ+1), ydummy

   ! Open file (note, strictly you should rewind when you open it)
   open (3, FILE='ptdata.txt', position='rewind', stat=ierror)
   if (ierror /= 0) then
      print *, 'Cannot open file'
      stop
   end if

   ! Read values
   do i = 1, 8192
      read (3, *, stat=ierror) ydummy, y(i)  
      if (ierror /= 0) then
         print *, 'Error or End of file found'
         stop
      end if
   end do

   ! Close file
   close (3)

   ! Do averaging.
   do i=1, size(y)-WINSIZ+1
      ! i is always index to first value in the window
      ! Next line does the averaging
      z(i) = sum(y(i:i+WINSIZ-1)) / real(WINSIZ)

      ! You can also do the averaging like this, but I would encourage
      ! you to become familiar with sum as it is more succint (and
      ! possibly quicker).
      ! z(i) = 0.0
      ! do j=i, i+WINSIZ-1
      !    z(i) = z(i) + y(j)
      ! end do
      ! z(i) = z(i) / real(WINSIZ)
   end do

   ! Print results
   do i=1, size(y)-WINSIZ+1
      print *, z(i)
   end do

end program

Thinking about it some more, the averaging can be made more efficient, since at each step all you do is remove the effect of the first point and add the effect of the next new point. So the averaging loop can be written as follows. The one subtraction and one addition at each step is nearly 4 times faster than the 8 additions in the loop above.

   ! Do averaging.
   z(1) = sum(y(1:WINSIZ))
   do i=2, size(y)-WINSIZ+1
      z(i) = z(i) - y(i-1) + y(i+WINSIZ-1)
   end do
   ! Convert sums to averages
   z = z / real(WINSIZ)
8 Jul 2011 11:55 #8535

Cheers! mate!

I finally worked out in the following way..

integer, parameter, dimension(*) :: a = [ 1, 4, 5, ..., 2 ]
integer                          :: i
real, dimension(10)              :: avg

avg = [ (sum(a(i * 4 + 1 : (i + 1) * 4)) / 4., i = 0, 9) ]
print *, avg
end
15 Jul 2011 12:43 #8578

The following variant of your code works with FTN95

real, dimension(40)    :: a 
integer                :: i 
real, dimension(10)    :: avg 
real*8 random
external random
!
do i = 1, size(a)
   a(i) = random()
end do
!
do i = 1,10
   avg(i) = sum (a(i*4-3 : i*4)) / 4.
end do
print *, avg 
!
avg = (/ (sum(a(i*4-3 : i*4)) / 4., i = 1, 10) /)
print *, avg 
end
Please login to reply.