SUBROUTINE filter_and_downsample(output, samples, downsample, filter, n, m)

  IMPLICIT NONE
  ! Variabili
  INTEGER m, n
  DOUBLE PRECISION, DIMENSION(n) :: samples
  DOUBLE PRECISION, DIMENSION(m) :: filter

  INTEGER downsample
  INTEGER k, i, s
  DOUBLE PRECISION, DIMENSION(n/downsample) :: output

!F2PY INTENT(IN) :: samples
!F2PY INTENT(IN) :: filter
!F2PY INTENT(IN) :: downsample
!F2PY INTENT(HIDE) :: n
!F2PY INTENT(HIDE) :: m
!F2PY INTENT(OUT) :: output

  ! Cominciamo a filtrare
  s = 1

  ! Applichiamo il filtro al sample k-esimo
  DO k = 1,n,downsample
     output(s) = 0
     i = 0
     DO WHILE( i < m .and. i < k )
        output(s) = output(s) + filter(i + 1)*samples(k - i)
        i = i + 1
     END DO

     ! Passiamo ai prossimi sample
     s = s + 1

  END DO

END SUBROUTINE

SUBROUTINE upsample_and_filter(output, samples, upsample, filter, n, m)

  ! dichiarazioni
  INTEGER :: m,n, upsample, i,j,s
  DOUBLE PRECISION, DIMENSION(n) :: samples
  DOUBLE PRECISION, DIMENSION(m) :: filter
  DOUBLE PRECISION, DIMENSION(upsample * n) :: output

  !F2PY INTENT(IN) samples
  !F2PY INTENT(IN) filter
  !F2PY INTENT(IN) upsample
  !F2PY INTENT(HIDE) m
  !F2PY INTENT(HIDE) n
  !F2PY INTENT(OUT) output

  s = 1
  DO i = 1, n
     ! Calcolo l'elemento in posizione s e s+1
     output(s) = 0
     output(s+1) = 0
     j = 0

     ! In questo ciclo calcoliamo la convoluzione del filtro 
     ! con il vettore upsampled sfruttando l'informazione che
     ! nei posti dispari c'è solo 0.
     DO WHILE(j < m .and. j < s)
        output(s)     = output(s) + samples(i - j/2) * filter(j + 1)
        output(s + 1) = output(s + 1) + samples(i - j/2) * filter(j + 2)
        j = j + 2
     END DO
     s = s + 2
  END DO

END SUBROUTINE
