program main use kind_defs use FourierTransform implicit none integer, parameter :: N = 2**15 ! 32768 complex(dbl) :: c(1:N) real(dbl) :: P(N), t(1:N), x(1:N), rand(1:N), x2(1:N), Ttot, f, limit integer :: i, ios, j ! Read the Data from the file: x = 0.0 open(50, file='Data/simpel_short.dat', status='old', action='read') do i = 1,N read(50, *, iostat=ios) t(i), x(i) if (ios/=0) exit enddo close(50) ! Add random noise to the data: call random_number(rand) x = x + 4*(rand-0.5) print *, "************************************************" print *, " Calculating spectrum of 'Data/simpel_short.dat'" print *, " plus random noise using FFT." print *, " Number of datapoints: ", N ! Calculate the FFT: c = fft(x) ! From this find the PowerSpectrum: P = abs(c)*sqrt(real(N)) ! Normering Som Matlab P = P/(real(N)/real(2)) ! Fysisk Normering ! Write the spectrum to file: Ttot = maxval(t)-minval(t) open(51, file='spectrum.dat', status='replace', action='write') do i = 1,N/2 f = (i-1)*(1e6/Ttot) ! The frequency in microHz write(51, '(2f20.8)') f, P(i) enddo close(51) print *, " Spectrum written to 'spectrum.dat'" print *, "************************************************" limit = 0.5 print *, " Removing all contributions with amplitude < ", limit print *, " This is done by setting the fouriertransform to 0" print *, " and doing the inverse FFT." ! Remove contributions from noise: do i = 1,N if (P(i) < limit) then c(i) = (0.0,0.0) endif enddo ! Transform back to filtered data: x2 = real(ifft(c)) ! Write the filtered data to file: open(51, file='filtered.dat', status='replace', action='write') do i = 1,N write(51, '(3f20.8)') t(i), x(i), x2(i) enddo close(51) print *, " Filtered data written to 'filtered.dat'." print *, "************************************************" end program