Line data    Source code

       1             : /**
       2             :  * \file ConvolutionFilter.cpp
       3             :  */
       4             : 
       5             : #include "ConvolutionFilter.h"
       6             : 
       7             : #include <cassert>
       8             : #include <algorithm>
       9             : 
      10             : namespace ATK
      11             : {
      12             :   template <typename DataType>
      13           1 :   ConvolutionFilter<DataType>::ConvolutionFilter()
      14           1 :     :Parent(1, 1)
      15             :   {
      16           1 :   }
      17             : 
      18             :   template<typename DataType_>
      19           1 :   void ConvolutionFilter<DataType_>::set_impulse(AlignedScalarVector impulse)
      20             :   {
      21           1 :     this->impulse = std::move(impulse);
      22           1 :     setup();
      23           1 :   }
      24             :   
      25             :   template<typename DataType_>
      26           1 :   void ConvolutionFilter<DataType_>::set_split_size(unsigned int split_size)
      27             :   {
      28           1 :     this->split_size = split_size;
      29           1 :     setup();
      30           1 :   }
      31             :   
      32             :   template<typename DataType_>
      33           3 :   void ConvolutionFilter<DataType_>::setup()
      34             :   {
      35           3 :     if(impulse.size() == 0 || split_size == 0)
      36             :     {
      37           2 :       return;
      38             :     }
      39             : 
      40           1 :     auto nb_splits = (impulse.size() + split_size - 1) / split_size;
      41           1 :     partial_frequency_input.assign(nb_splits - 1, AlignedComplexVector(split_size * 2, 0));
      42             :     // Pad with zeros so the convolution is easier created.
      43           1 :     impulse.resize((partial_frequency_input.size() + 1) * split_size, 0);
      44             : 
      45           1 :     temp_out_buffer.assign(split_size * 2, 0);
      46           1 :     processor.set_size(split_size * 2);
      47             :     
      48             :     // The size is twice as big than the impulse, less
      49           1 :     partial_frequency_impulse.assign(split_size * 2 * (nb_splits - 1), 0);
      50          79 :     for(gsl::index i = 1; i < nb_splits; ++i)
      51             :     {
      52          78 :       processor.process_forward(impulse.data() + i * split_size, partial_frequency_impulse.data() + (i - 1) * split_size * 2, split_size);
      53             :     }
      54             : 
      55           1 :     input_delay = split_size - 1;
      56             :     
      57           1 :     Parent::setup();
      58             :   }
      59             : 
      60             :   template<typename DataType_>
      61       16376 :   void ConvolutionFilter<DataType_>::compute_convolutions() const
      62             :   {
      63     2112500 :     for(gsl::index i = 0; i < split_size; ++i)
      64             :     {
      65     2096130 :       temp_out_buffer[i] = temp_out_buffer[i + split_size];
      66     2096130 :       temp_out_buffer[i + split_size] = 0;
      67             :     }
      68             :     
      69       16376 :     result.assign(2 * split_size, 0);
      70       16376 :     DataType* ATK_RESTRICT result_ptr_orig = reinterpret_cast<DataType*>(result.data());
      71       16376 :     const DataType* ATK_RESTRICT partial_frequency_impulse_ptr_orig = reinterpret_cast<const DataType*>(partial_frequency_impulse.data());
      72             : 
      73             :     // offset in the impulse frequencies
      74       16376 :     int64_t offset = 0;
      75     1293700 :     for(const auto& buffer: partial_frequency_input)
      76             :     {
      77     1277330 :       DataType* ATK_RESTRICT result_ptr = result_ptr_orig;
      78     1277330 :       const DataType* ATK_RESTRICT buffer_ptr = reinterpret_cast<const DataType*>(buffer.data());
      79     1277330 :       const DataType* ATK_RESTRICT partial_frequency_impulse_ptr = partial_frequency_impulse_ptr_orig + offset;
      80             :       // Add the frequency result of this partial FFT
      81   164775000 :       for(gsl::index i = 0; i < split_size; ++i)
      82             :       {
      83   163498000 :         DataType br1 = *(buffer_ptr++);
      84   163498000 :         DataType bi1 = *(buffer_ptr++);
      85   163498000 :         DataType pr1 = *(partial_frequency_impulse_ptr++);
      86   163498000 :         DataType pi1 = *(partial_frequency_impulse_ptr++);
      87   163498000 :         DataType br2 = *(buffer_ptr++);
      88   163498000 :         DataType bi2 = *(buffer_ptr++);
      89   163498000 :         DataType pr2 = *(partial_frequency_impulse_ptr++);
      90   163498000 :         DataType pi2 = *(partial_frequency_impulse_ptr++);
      91   163498000 :         *(result_ptr++) += br1*pr1-bi1*pi1;
      92   163498000 :         *(result_ptr++) += br1*pi1+pr1*bi1;
      93   163498000 :         *(result_ptr++) += br2*pr2-bi2*pi2;
      94   163498000 :         *(result_ptr++) += br2*pi2+pr2*bi2;
      95             :       }
      96     1277330 :       offset += 4 * split_size;
      97             :     }
      98             : 
      99       32752 :     std::vector<DataType> ifft_result(2*split_size, 0);
     100       16376 :     processor.process_backward(result.data(), ifft_result.data(), 2*split_size);
     101     4208630 :     for(gsl::index i = 0; i < 2*split_size; ++i)
     102             :     {
     103     4192260 :       temp_out_buffer[i] += ifft_result[i] * split_size * 2;
     104             :     }
     105       16376 :   }
     106             : 
     107             :   template<typename DataType_>
     108       16376 :   void ConvolutionFilter<DataType_>::process_new_chunk(int64_t position) const
     109             :   {
     110       16376 :     if(partial_frequency_input.empty())
     111             :     {
     112           0 :       return;
     113             :     }
     114       16376 :     partial_frequency_input.pop_back();
     115       32752 :     AlignedComplexVector chunk(2 * split_size);
     116       16376 :     processor.process_forward(converted_inputs[0] + position - split_size, chunk.data(), split_size);
     117             :     
     118       16376 :     partial_frequency_input.push_front(std::move(chunk));
     119       16376 :     compute_convolutions();
     120             :   }
     121             : 
     122             :   template<typename DataType_>
     123       18408 :   void ConvolutionFilter<DataType_>::process_impulse_beginning(int64_t processed_size, unsigned int size_to_process) const
     124             :   {
     125       18408 :     const DataType* ATK_RESTRICT input = converted_inputs[0] + processed_size;
     126       18408 :     const DataType* ATK_RESTRICT impulse_ptr = impulse.data();
     127       18408 :     DataType* ATK_RESTRICT output = outputs[0] + processed_size;
     128             : 
     129     2114540 :     for(gsl::index i = 0; i < size_to_process; ++i)
     130             :     {
     131     2096130 :       output[i] = temp_out_buffer[split_position + i];
     132             :     }
     133             : 
     134     2374630 :     for(gsl::index j = 0; j < input_delay + 1; ++j)
     135             :     {
     136   270661000 :       for(gsl::index i = 0; i < size_to_process; ++i)
     137             :       {
     138   268304000 :         output[i] += impulse_ptr[j] * input[static_cast<int>(i) - static_cast<int>(j)];
     139             :       }
     140             :     }
     141       18408 :   }
     142             : 
     143             :   template<typename DataType_>
     144        2047 :   void ConvolutionFilter<DataType_>::process_impl(gsl::index size) const
     145             :   {
     146        2047 :     assert(input_sampling_rate == output_sampling_rate);
     147             :     
     148        2047 :     unsigned int processed_size = 0;
     149       16361 :     do
     150             :     {
     151             :       // We can only process split_size elements at a time, but if we already have some elements in the buffer,
     152             :       // we need to take them into account.
     153       18408 :       unsigned int size_to_process = std::min(split_size - split_position, static_cast<unsigned int>(size) - processed_size);
     154             : 
     155       18408 :       process_impulse_beginning(processed_size, size_to_process);
     156             : 
     157       18408 :       split_position += size_to_process;
     158       18408 :       processed_size += size_to_process;
     159       18408 :       if(split_position == split_size)
     160             :       {
     161       16376 :         process_new_chunk(processed_size);
     162       16376 :         split_position = 0;
     163             :       }
     164       18408 :     }while(processed_size != size);
     165        2047 :   }
     166             : 
     167             :   template class ConvolutionFilter<double>;
     168             : }