compressed_sensing.F90

Go to the documentation of this file.

!! Copyright (C) 2011 X. Andrade
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
 
#include <global.h>
 
module compressed_sensing_oct_m
  use blas_oct_m
  use bpdn_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
 
  ! these values are copied from td/spectrum.f90
  integer, parameter ::            &
    SPECTRUM_TRANSFORM_LAPLACE = 1,  &
    spectrum_transform_sin     = 2,  &
    spectrum_transform_cos     = 3
 
  public ::                                        &
    compressed_sensing_t,                          &
    compressed_sensing_init,                       &
    compressed_sensing_end,                        &
    compressed_sensing_spectral_analysis
 
  type compressed_sensing_t
    private
    real(real64)   :: sigma
    integer :: ntime
    real(real64)   :: dtime
    real(real64)   :: stime
    integer :: nfreq
    real(real64)   :: dfreq
    real(real64)   :: sfreq
    type(bpdn_matrix) :: fourier_matrix
  end type compressed_sensing_t
 
contains
 
  subroutine compressed_sensing_init(this, transform_type, ntime, dtime, stime, nfreq, dfreq, sfreq, noise)
    type(compressed_sensing_t),  intent(out) :: this
    integer,                     intent(in)  :: transform_type
    integer,                     intent(in)  :: ntime
    real(real64),                intent(in)  :: dtime
    real(real64),                intent(in)  :: stime
    integer,                     intent(in)  :: nfreq
    real(real64),                intent(in)  :: dfreq
    real(real64),                intent(in)  :: sfreq
    real(real64),                intent(in)  :: noise
 
    integer :: itime, ifreq, type
    real(real64)   :: time, freq
 
    push_sub(compressed_sensing_init)
 
    this%sigma = noise
 
    this%ntime = ntime
    this%dtime = dtime
    this%stime = stime
    this%nfreq = nfreq
    this%dfreq = dfreq
    this%sfreq = sfreq
 
    if (transform_type == spectrum_transform_laplace) then
 
      call bpdn_matrix_init(this%fourier_matrix, this%ntime, this%nfreq, explicit_matrix)
 
      do ifreq = 1, this%nfreq
        freq = (ifreq - 1)*this%dfreq + this%sfreq
 
        do itime = 1, this%ntime
          time = (itime - 1)*this%dtime + this%stime
          this%fourier_matrix%matrix(itime, ifreq) = exp(-freq*time)
        end do
 
      end do
 
    else
 
      select case (transform_type)
      case (spectrum_transform_sin)
        type = sin_matrix
      case (spectrum_transform_cos)
        type = cos_matrix
      end select
 
      call bpdn_matrix_init(this%fourier_matrix, this%ntime, this%nfreq, type)
      call bpdn_matrix_set_delta(this%fourier_matrix, this%dtime, this%dfreq)
 
    end if
 
    pop_sub(compressed_sensing_init)
  end subroutine compressed_sensing_init
 
  ! -------------------------------------------------------------------
 
  subroutine compressed_sensing_end(this)
    type(compressed_sensing_t),  intent(inout) :: this
 
    push_sub(compressed_sensing_end)
 
    call bpdn_matrix_end(this%fourier_matrix)
 
    pop_sub(compressed_sensing_end)
  end subroutine compressed_sensing_end
 
  ! -------------------------------------------------------------------
 
  subroutine compressed_sensing_spectral_analysis(this, time_function, freq_function)
    type(compressed_sensing_t),  intent(inout) :: this
    real(real64),                intent(in)    :: time_function(:)
    real(real64),                intent(out)   :: freq_function(:)
 
    integer :: ierr
    real(real64), allocatable :: tf_normalized(:)
    real(real64) :: nrm
 
    push_sub(compressed_sensing_spectral_analysis)
 
    safe_allocate(tf_normalized(1:this%ntime))
 
    ! to avoid numerical problems we work with a normalized rhs
    nrm = dnrm2(this%ntime, time_function(1), 1)
 
    if (nrm > 1e-8_real64) then
      tf_normalized(1:this%ntime) = time_function(1:this%ntime)/nrm
    else
      tf_normalized(1:this%ntime) = time_function(1:this%ntime)
      nrm = m_one
    end if
 
    call bpdn(this%ntime, this%nfreq, this%fourier_matrix, tf_normalized, this%sigma, freq_function, ierr, activesetit = 50)
 
    safe_deallocate_a(tf_normalized)
 
    ! scale by the missing factors
    freq_function(1:this%nfreq) = nrm/(this%dfreq*m_two/m_pi)*freq_function(1:this%nfreq)
 
    if (ierr < 0) then
      message(1) = 'The Basis Pursuit Denoising process failed to converge.'
      call messages_warning(1)
    end if
 
    pop_sub(compressed_sensing_spectral_analysis)
  end subroutine compressed_sensing_spectral_analysis
 
end module compressed_sensing_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: