Chris@1267: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@1267: 
Chris@1267: /*
Chris@1267:     Sonic Visualiser
Chris@1267:     An audio file viewer and annotation editor.
Chris@1267:     Centre for Digital Music, Queen Mary, University of London.
Chris@1267:     This file copyright 2006-2016 Chris Cannam and QMUL.
Chris@1267:     
Chris@1267:     This program is free software; you can redistribute it and/or
Chris@1267:     modify it under the terms of the GNU General Public License as
Chris@1267:     published by the Free Software Foundation; either version 2 of the
Chris@1267:     License, or (at your option) any later version.  See the file
Chris@1267:     COPYING included with this distribution for more information.
Chris@1267: */
Chris@1267: 
Chris@1267: #include "ColumnOp.h"
Chris@1267: 
Chris@1267: #include <cmath>
Chris@1267: #include <algorithm>
Chris@1267: #include <iostream>
Chris@1267: 
Chris@1283: #include "base/Debug.h"
Chris@1283: 
Chris@1267: using namespace std;
Chris@1267: 
Chris@1267: ColumnOp::Column
Chris@1267: ColumnOp::fftScale(const Column &in, int fftSize)
Chris@1267: {
Chris@1267:     return applyGain(in, 2.0 / fftSize);
Chris@1267: }
Chris@1267: 
Chris@1267: ColumnOp::Column
Chris@1267: ColumnOp::peakPick(const Column &in)
Chris@1267: {
Chris@1837:     Column out(in.size(), 0.f);
Chris@1267: 
Chris@1267:     for (int i = 0; in_range_for(in, i); ++i) {
Chris@1267:         if (isPeak(in, i)) {
Chris@1267:             out[i] = in[i];
Chris@1267:         }
Chris@1267:     }
Chris@1267:     
Chris@1267:     return out;
Chris@1267: }
Chris@1267: 
Chris@1267: ColumnOp::Column
Chris@1267: ColumnOp::normalize(const Column &in, ColumnNormalization n) {
Chris@1267: 
Chris@1267:     if (n == ColumnNormalization::None || in.empty()) {
Chris@1267:         return in;
Chris@1267:     }
Chris@1394:     
Chris@1394:     float shift = 0.f;
Chris@1394:     float scale = 1.f;
Chris@1267: 
Chris@1394:     if (n == ColumnNormalization::Range01) {
Chris@1394: 
Chris@1394:         float min = 0.f;
Chris@1394:         float max = 0.f;
Chris@1394:         bool have = false;
Chris@1394:         for (auto v: in) {
Chris@1394:             if (v < min || !have) {
Chris@1394:                 min = v;
Chris@1394:             }
Chris@1394:             if (v > max || !have) {
Chris@1394:                 max = v;
Chris@1394:             }
Chris@1394:             have = true;
Chris@1394:         }
Chris@1394:         if (min != 0.f) {
Chris@1394:             shift = -min;
Chris@1394:             max -= min;
Chris@1394:         }
Chris@1394:         if (max != 0.f) {
Chris@1394:             scale = 1.f / max;
Chris@1394:         }
Chris@1394: 
Chris@1394:     } else if (n == ColumnNormalization::Sum1) {
Chris@1267: 
Chris@1267:         float sum = 0.f;
Chris@1267: 
Chris@1267:         for (auto v: in) {
Chris@1267:             sum += fabsf(v);
Chris@1267:         }
Chris@1267: 
Chris@1267:         if (sum != 0.f) {
Chris@1267:             scale = 1.f / sum;
Chris@1267:         }
Chris@1267: 
Chris@1267:     } else {
Chris@1267: 
Chris@1267:         float max = 0.f;
Chris@1267: 
Chris@1267:         for (auto v: in) {
Chris@1267:             v = fabsf(v);
Chris@1267:             if (v > max) {
Chris@1267:                 max = v;
Chris@1267:             }
Chris@1267:         }
Chris@1267: 
Chris@1267:         if (n == ColumnNormalization::Max1) {
Chris@1267:             if (max != 0.f) {
Chris@1267:                 scale = 1.f / max;
Chris@1267:             }
Chris@1267:         } else if (n == ColumnNormalization::Hybrid) {
Chris@1267:             if (max > 0.f) {
Chris@1267:                 scale = log10f(max + 1.f) / max;
Chris@1267:             }
Chris@1267:         }
Chris@1267:     }
Chris@1267: 
Chris@1394:     return applyGain(applyShift(in, shift), scale);
Chris@1267: }
Chris@1267: 
Chris@1267: ColumnOp::Column
Chris@1267: ColumnOp::distribute(const Column &in,
Chris@1267:                      int h,
Chris@1267:                      const vector<double> &binfory,
Chris@1267:                      int minbin,
Chris@1267:                      bool interpolate)
Chris@1267: {
Chris@1837:     Column out(h, 0.f);
Chris@1267:     int bins = int(in.size());
Chris@1267: 
Chris@1303:     if (interpolate) {
Chris@1303:         // If the bins are all closer together than the target y
Chris@1303:         // coordinate increments, then we don't want to interpolate
Chris@1303:         // after all. But because the binfory mapping isn't
Chris@1303:         // necessarily linear, just checking e.g. whether bins > h is
Chris@1303:         // not enough -- the bins could still be spaced more widely at
Chris@1303:         // either end of the scale. We are prepared to assume however
Chris@1303:         // that if the bins are closer at both ends of the scale, they
Chris@1303:         // aren't going to diverge mysteriously in the middle.
Chris@1303:         if (h > 1 &&
Chris@1303:             fabs(binfory[1] - binfory[0]) >= 1.0 &&
Chris@1303:             fabs(binfory[h-1] - binfory[h-2]) >= 1.0) {
Chris@1303:             interpolate = false;
Chris@1303:         }
Chris@1303:     }
Chris@1303:     
Chris@1267:     for (int y = 0; y < h; ++y) {
Chris@1267: 
cannam@1301:         if (interpolate) {
Chris@1267: 
Chris@1267:             double sy = binfory[y] - minbin - 0.5;
Chris@1267:             double syf = floor(sy);
Chris@1267: 
Chris@1267:             int mainbin = int(syf);
Chris@1267:             int other = mainbin;
Chris@1267:             if (sy > syf) {
Chris@1267:                 other = mainbin + 1;
Chris@1267:             } else if (sy < syf) {
Chris@1267:                 other = mainbin - 1;
Chris@1267:             }
Chris@1267: 
Chris@1267:             if (mainbin < 0) {
Chris@1267:                 mainbin = 0;
Chris@1267:             }
Chris@1267:             if (mainbin >= bins) {
Chris@1267:                 mainbin = bins - 1;
Chris@1267:             }
Chris@1267: 
Chris@1267:             if (other < 0) {
Chris@1267:                 other = 0;
Chris@1267:             }
Chris@1267:             if (other >= bins) {
Chris@1267:                 other = bins - 1;
Chris@1267:             }
Chris@1267: 
Chris@1267:             double prop = 1.0 - fabs(sy - syf);
Chris@1267:             
Chris@1267:             double v0 = in[mainbin];
Chris@1267:             double v1 = in[other];
Chris@1267:                 
Chris@1267:             out[y] = float(prop * v0 + (1.0 - prop) * v1);
Chris@1267: 
Chris@1267:         } else {
Chris@1267:             
Chris@1267:             double sy0 = binfory[y] - minbin;
Chris@1267: 
Chris@1267:             double sy1;
Chris@1267:             if (y+1 < h) {
Chris@1267:                 sy1 = binfory[y+1] - minbin;
Chris@1267:             } else {
Chris@1267:                 sy1 = bins;
Chris@1267:             }
Chris@1267: 
Chris@1267:             int by0 = int(sy0 + 0.0001);
Chris@1267:             int by1 = int(sy1 + 0.0001);
Chris@1283: 
Chris@1283:             if (by0 < 0 || by0 >= bins || by1 > bins) {
Chris@1283:                 SVCERR << "ERROR: bin index out of range in ColumnOp::distribute: by0 = " << by0 << ", by1 = " << by1 << ", sy0 = " << sy0 << ", sy1 = " << sy1 << ", y = " << y << ", binfory[y] = " << binfory[y] << ", minbin = " << minbin << ", bins = " << bins << endl;
Chris@1283:                 continue;
Chris@1283:             }
Chris@1267:                 
Chris@1267:             for (int bin = by0; bin == by0 || bin < by1; ++bin) {
Chris@1267: 
Chris@1267:                 float value = in[bin];
Chris@1267: 
Chris@1267:                 if (bin == by0 || value > out[y]) {
Chris@1267:                     out[y] = value;
Chris@1267:                 }
Chris@1267:             }
Chris@1267:         }
Chris@1267:     }
Chris@1267: 
Chris@1267:     return out;
Chris@1267: }