wegen/sauSynKlick/SynKlick.cc

/*
 SynKlick  -  repair synchronous skips in raw sound recording

 (C) Ichthyostega (Hermann Vosseler), 2014

All code in this package is provided "as is" and WITHOUT ANY WARRANTY.
It may be copied and used under the terms and conditions of the
GNU PUBLIC LICENCE (GPL) Version 2 and above.

*********************************************************************/

#include <sndfile.h>
#include <stdio.h>

#include <string>
#include <vector>
#include <iostream>
#include <cassert>
#include <cstdlib>
#include <cmath>

using std::string;
using std::vector;
using std::cout;
using std::cerr;
using std::max;
using std::abs;

const size_t CHUNK_SIZE  = 20000; //44100;  ///TODO for debugging
const size_t CONTEXT    = 100;
const size_t BUFF_SIZE = CHUNK_SIZE + 2*CONTEXT;

const float EPSILON = 1.0f / (1 << 23);      // one normalised sample step for 24bit samples 
const float SILENCE_GATE = 4096.0f * EPSILON;// signal below this level (-33dB) is too low to detect jumps
const float TRIGGER_FACTOR = 5.0f;           // discontinuity detected by local change compared to average change
const float CONFIRMATION_DEGREE = 0.66f;     // degree of confirmation by multiple tracks for an event to be accepted
const uint  AVERAGE_SAMPLES = 50;            // number of past samples to use for average delta

const bool VERBOSE = true;

/* ==== Helpers ==== */

void
error (string msg)
{
  cerr << "processing FAILED: "
       << msg << "\n\n";
  exit(EXIT_FAILURE);
}

string
name_without_path (string rawName)
{
  auto pos = rawName.rfind('/');
  return (pos == string::npos)? rawName
                              : rawName.substr(pos+1);
}

void
reduceBufferPos (size_t & pos, size_t elementsToRemove)
{
  if (pos < elementsToRemove)
    pos = 0;
  else
    pos -= elementsToRemove;
}




/* ==== Abstraction: single channel to be processed ==== */

/**
 * Helper to decide if a glitch event should be considered,
 * while visiting multiple tracks at the same location.
 */
class EventScore
  {
    uint contributions_ = 0;
    uint positives_ = 0;
    float prevDeg_ = 0.0f;
    
  public:
    float
    confirmationDegree()
      {
        return contributions_? float(positives_)/contributions_
                             : 0.0f;
      }
    
    operator bool()
      {
        return CONFIRMATION_DEGREE <= this->confirmationDegree();
      }
    
    void
    recordConfirmation(bool confirmed)
      {
        prevDeg_ = confirmationDegree(); // record for diagnostics
        
        ++contributions_;
        if (confirmed)
          ++positives_;
      }
    
    string
    showLastUpdate()
      {
        float deg = confirmationDegree();
        return deg == prevDeg_? "neutral"
                              : (deg > prevDeg_? "conf "
                                               : "rej  "
                                ) + std::to_string(deg);
      }
  };



class Channel
  {
    SNDFILE *in_ = nullptr,
            *out_ = nullptr;
    SF_INFO info_in_,
            info_out_;
    
    float buff_[BUFF_SIZE];
    size_t numSamples_ = 0;
    size_t workPos_    = 0;
    float ema_ = 1.0f / TRIGGER_FACTOR;
    
    
  public:
    string name_;
    static bool hitEnd;
    
    void
    open (string name, string target)
      {
        this->name_ = name_without_path(name);
        target += name_;
        
        info_in_.format = 0;
        in_ = sf_open (name.c_str(), SFM_READ, &info_in_);
        if (!in_)
          error ("could not open "+name+" -- " + sf_strerror (NULL));
        if (1 != info_in_.channels)
          error ("only Mono streams supported -- but '"+name_
                +"' was found to have "+std::to_string(info_in_.channels)+" channels.");
        
        info_out_.format = info_in_.format;
        info_out_.channels = info_in_.channels;
        info_out_.samplerate = info_in_.samplerate;
        out_ = sf_open(target.c_str(), SFM_WRITE, &info_out_);
        if (!out_)
          error ("unable to create target "+target+" -- " + sf_strerror (NULL));
      }
    
    
   ~Channel()
      {
        if (in_) sf_close (in_);
        if (out_)
          if (0 != sf_close (out_))
            cerr << "WARNING: problems while closing target '" << name_
                 << "' -- " << sf_strerror(out_) << "\n";
        
        in_ = out_ = nullptr;
      }
    
    
    void
    fillBuffer()
      {
        assert(numSamples_ < BUFF_SIZE);
        float* startPos = buff_ + numSamples_;
        size_t samplesToRead = BUFF_SIZE - numSamples_;
        numSamples_ = sf_read_float(in_, startPos, samplesToRead);
        if (numSamples_ < samplesToRead)
          hitEnd = true; // causes the processing to end with this chunk.
        
        numSamples_ += (startPos - buff_);
        assert(workPos_ <= numSamples_);
      }
    
    void
    writeBuffer()
      {
        size_t samplesToWrite = numSamples_ - (hitEnd? 0 : 2*CONTEXT);
        size_t actuallyWritten = sf_write_float(out_, buff_, samplesToWrite);
        if (actuallyWritten < samplesToWrite)
          cerr << "Alarm: only "<<actuallyWritten<<" instead of "<<samplesToWrite
               <<" samples written for " << name_ <<"\n";
        
        if (!hitEnd)
          // move processing context samples to the beginning of the buffer;
          // they will be written together with the next chunk of data
          {
            for (uint i=0; i < 2*CONTEXT; ++i)
              buff_[i] = buff_[i + actuallyWritten];
            reduceBufferPos(numSamples_, actuallyWritten);
            reduceBufferPos(workPos_,    actuallyWritten);
          }
        assert(workPos_ <= numSamples_);
      }
    
    /**
     * investigate a single sample position in the current buffer content.
     * The delta to the preceding sample(s) is used to detect a signal glitch.
     * @param i position of current sample
     * @param event state of other tracks investigated thus far
     * @return the event detection state including the results of investigating this sample.
     *         If the sample is positively an glitch event, and the state is undecided thus far,
     *         the event state is confirmed; if a glitch can not be confirmed, the state is
     *         cleared. Otherwise, the undecided state is passed on.
     */
    void
    visit (uint i, EventScore& event)
      {
        if (i + CONTEXT < numSamples_
            && CONTEXT  <= i
            && workPos_ <= i)
          {
            workPos_ = i;
            float averageDelta = updateMovingAverage();
            float changeFactor = delta() / (averageDelta + EPSILON);
            if (SILENCE_GATE < sampleLevel())
              event.recordConfirmation (changeFactor > TRIGGER_FACTOR);
          }
      }
    
  private:
    float
    sampleLevel()
      {
        auto curr = abs(buff_[workPos_]);
        auto prev = abs(buff_[workPos_ - 1]);
        
        return 0 < workPos_? max(prev, curr)
                           : curr;
      }
    
    float
    delta()
      {
        if (0 < workPos_)
          return abs(buff_[workPos_] - buff_[workPos_-1]);
        else
          return 0;
      }
    
    float
    updateMovingAverage()
      {
        return ema_ = (delta() + ema_*(AVERAGE_SAMPLES - 1)) / AVERAGE_SAMPLES;
      }
    
    float*
    accessCurrentContext()
      {
        assert(CONTEXT <= workPos_);
        assert(CONTEXT + workPos_ <= numSamples_);
        return buff_ + (workPos_ - CONTEXT);
      }
    friend class Patch;
    
  };



/**
 * Fix a glitch detected in the corresponding channel.
 * When the Patch is applied, internal work buffers are used
 * to create two (partially synthetic) copies of the channel's signal
 * which are then cross faded. The result is pushed back into the
 * output of the underlying channel.
 */
class Patch
  {
    static const size_t LEN = 2 * CONTEXT;
    
    float workBuff[LEN];
    float* rawBuff;
    
    
    void
    reconstructMirroredSignal()
      {
        for (uint i=0; i < LEN; ++i)
          workBuff[i] = rawBuff[LEN-1 - i];
      }
    
    
  public:
    Patch (Channel& c)
      : rawBuff{c.accessCurrentContext()}
      {
        reconstructMirroredSignal();
      }
    
    void
    apply()
      {
        for (uint i=0; i < LEN; ++i)
          {
            float fade = (i<CONTEXT) ? CONTEXT - (i+1)
                                     : i - CONTEXT;
            fade = 1.0 - fade/CONTEXT;
            // i=-1        -> 0
            // i=CONTEXT-1 -> 1
            // i=CONTEXT   -> 1
            // i=2*CONTEXT -> 0
            
            float p = fade * 0.5f,
                  q = 1.0 - p;
            
            rawBuff[i] = q*rawBuff[i] + p*workBuff[i];
          }
      }
  };





/* ==== high-level operation ==== */

int num_sources = 0;
vector<Channel> streams(0);
bool Channel::hitEnd = false;


void
openStreams (char* name[], string targetDir)
{
  assert (num_sources);
  assert (0 == streams.size());
  streams.resize(num_sources);
  
  assert (targetDir.length());
  if (targetDir.back() != '/')
    targetDir += '/';
  
  for (uint i=0; i < num_sources; ++i)
    streams[i].open (name[i], targetDir);
}

void
closeStreams()
{
  streams.clear();
}


void
eventConfirmationDiagnostic(EventScore& event, string name, size_t absPos)
{
    static size_t trail = 0;
    if (trail && trail != absPos) trail=0;
    if (event && !trail) trail = absPos;
    
    if (trail)
      cerr <<absPos<<"\t jump? "
           <<name<<"\t :: "
           <<event.showLastUpdate()
           <<"\n";
}

void
process()
{
  size_t absStartPos = 0; // for diagnostic messages
  
  while (!Channel::hitEnd)
    {
      for (Channel& channel : streams)
        channel.fillBuffer();
      
      for (uint i=0; i < BUFF_SIZE; ++i)
        {
          EventScore event;
          for (Channel& channel : streams)
            {
              channel.visit(i, event);
              if (VERBOSE)
                eventConfirmationDiagnostic (event, channel.name_, absStartPos+i);
            }
          if (event)
            {
              cout << "PATCH at: "<< absStartPos+i <<"\n";
              for (Channel& channel : streams)
                Patch(channel).apply();
            }
        }
      
      for (Channel& channel : streams)
        channel.writeBuffer();
      
      absStartPos += CHUNK_SIZE;
    }
}





int
main (int argc, char* argv[]) {
    if (argc < 3 || string("-h") == argv[1])
      {
        cerr << "SynKlick <INPUTFILE> [<INPUTFILE2> ...] <TARGETDIR>\n"
             << "\nRepair a very special kind of defects in sound recordings:\n"
             << "synchronous skips in a set of input sound files. Possible events\n"
             << "will be detected by scanning the deltas between neighbouring samples\n"
             << "and confirmed by jumps at the same position in all files.\n"
             << "The jump will be amended by a short cross-blend of about 4.5ms.\n"
             << "\nA fixed copy of each file will be written to <TARGETDIR>\n\n";
        
        return EXIT_FAILURE;
      }
    num_sources = argc - 2;
    string targetDir = argv[argc-1];
    
    openStreams (&argv[1], targetDir);
    process();
    closeStreams();
    cout << "\nSUCCESS\n"
         << num_sources << " files read and processed --> "
         << targetDir << "\n\n";
    return EXIT_SUCCESS;
}