 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
Srabon
Joined: 22 Feb 2011
Posts: 14
|
 Posted: Mon Jul 04, 2011 10:05 pm Post subject: code to window averaging filter |
 |
|
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!
| Code: | 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] |
|
| Back to top |
|
 |
davidb
Joined: 17 Jul 2009
Posts: 560
Location: UK
|
| Posted: Thu Jul 07, 2011 10:56 pm Post subject: Re: code to window averaging filter |
|
|
| Srabon wrote: | 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 
I provide a more optimsed version below.
| Code: |
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.
| Code: |
! 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)
|
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl |
|
| Back to top |
|
|
Srabon
Joined: 22 Feb 2011
Posts: 14
|
| Posted: Sat Jul 09, 2011 12:55 am Post subject: |
|
|
Cheers! mate!
I finally worked out in the following way..
| Code: |
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 |
|
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2615
Location: Sydney
|
| Posted: Fri Jul 15, 2011 1:43 am Post subject: |
|
|
The following variant of your code works with FTN95
| Code: | 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 |
|
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
