package de.maramuse.soundcomp.effect;
/*
* Copyright 2010 Jan Schmidt-Reinisch
*
* SoundComp - a sound processing library
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; in version 2.1
*
* This library 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*
*/
/**
* An implementation of a frequency scaling by TDHS (time domain harmonic scaling)
*/
import de.maramuse.soundcomp.filter.StateVariableFilter;
import de.maramuse.soundcomp.process.ConstStream;
import de.maramuse.soundcomp.process.NamedSource;
import de.maramuse.soundcomp.process.ParameterMap;
import de.maramuse.soundcomp.process.ProcessElement;
import de.maramuse.soundcomp.process.SourceStore;
import de.maramuse.soundcomp.process.StandardParameters;
import de.maramuse.soundcomp.process.Stateful;
import de.maramuse.soundcomp.process.TypeMismatchException;
import de.maramuse.soundcomp.process.UnknownConnectionException;
import de.maramuse.soundcomp.process.ValueType;
import de.maramuse.soundcomp.process.StandardParameters.Parameter;
import de.maramuse.soundcomp.util.GlobalParameters;
import de.maramuse.soundcomp.util.ReadOnlyMap;
import de.maramuse.soundcomp.util.ReadOnlyMapImpl;
import de.maramuse.soundcomp.util.NativeObjects;
/**
* For better decision of which way to go, we will implement a few planned modes of operation:
*
* Mode 0:
*
* Don't care about artifacts. Just do the transformation and let the buffer over/underflows happen. This is the
* cheapest way and its base algorithm will be used in any following mode. The buffer handling here will be: Determine
* the resampling factor F (f_out/f_in). For a semitone upwards, this will e.g. be 2^(1/12). Reject any resampling
* factors outside (0.1 .. 2). Split this number F into 2 factors, one (U, fractional) for upsampling and one (D,
* natural) for downsampling, so that F=D/U and D = 10. In our example, D=10, U=10*2^(-1/12)
*
* For each D samples that we will finally emit, we must insert U samples into the pipeline. Since U is not natural, we
* must do some error distribution there: some times we have to insert INT(U) samples, sometimes INT(U+1). We will
* aggregate the frational error to decide which case is relevant for each input sample. This upsampled data stream is
* LP-filtered with a SVF with an edge frequency of SAMPLERATE/2/D (actually the filter runs on a data stream D times as
* fast, so this simulates filtering on SAMPLERATE/2 since the filter does not know about the faked higher sample rate).
* Then, from the resulting upsampled and filtered data stream, take every D-th sample as final output.
*
* For F<1, this leads to buffer underruns, since we fill the pipe with less data than we attempt to get from it. For
* F>1, the opposite will happen: the pipe content will grow. Since we use a FIFO ring buffer, both cases will equally
* lead to "jumps" in the input stream by the size of the ring buffer. Mode 0 ignores these glitches, they will be
* addressed in following modes.
*
* Mode 1:
*
* split the input signal into two paths. one is double band pass filtered via state variable filters at an edge
* frequency of 20 and 50 hz, and then used for edge detection (upward zero crossings) to get a period duration D (in
* samples) hint. This should be largely cheaper than the usual autocorrelation mechanisms, and in most cases at least
* as precise (unless we have periodic signals, where autocorrelation is brilliant. But periodic signals will be more an
* exception here).
*
* the other path is upsampled by a large scale factor S (with S/F being natural and s<=100) then, depending on whether
* to shift up or down, either S*D sized fragments are cut out with windowing, or they are duplicated with windowing,
* from a buffer size of S*D/abs(F-1). After this, the resulting signal is low pass filtered by half the original
* sampling frequency, then every S/F-th sample is taken out to form the output sample stream.
*
* The upsampling may be a by non-natural factor. It is simply done by repeating the input value the necessary number of
* times. In case the factor is non-natural, it is approximated by spreading the error between repeating the ceil and
* floor values alternatively. This is really no problem since the resulting error is mitigated by the anti alias
* filtering that comes afterwards. It is essential that the downsampling is carried out by a natural factor since no
* error spreading is possible with post processing, because the downsampling is the last step.
*
* We will not allow any base frequency below audible range (more precisely, not even below SAMPLINGRATE/100).
* Otherwise, the buffer lengths would get way too long, and the delays would get really annoying. Therefore, the edge
* detection will not look at larger portions than 100 samples. Should the edge detection fail, the largest allowed
* buffer size is applied.
*
* Note that fast changes in the frequency ratio will almost surely produce artefacts since the pointer arithmetics in
* the up/downsampling buffers are not prepared for this parameter to change.
*
*
*
* Mode 2:
*
* Do a regular covariance detection. Limit this covariance to a relatively short range, but search in a longer range.
* So we can find shorter similar curve fragments within a signal of longer period without having to do the expensive
* covariance detection over thousands of samples. A range of 30 samples of best covariance should suffice, but search
* this within a range of at least 2000 samples (SAMPLERATE/20Hz). This method is more expensive, but I expect it to
* better suppress artifacts. As mode 1, it does not necessarily require a periodic input. My favourite candidate.
*
* Mode 3:
*
*
* @author jssr67
*
*/
public class TDHS implements ProcessElement, Stateful {
/**
* MODE determines the algorithm that is used for finding suitable "cutting points" in the sample stream. See
* explanation above
*/
private static int MODE=1;
private static final ReadOnlyMapImpl<Integer, ValueType> recSourceMap=new ReadOnlyMapImpl<Integer, ValueType>();
private static final ReadOnlyMapImpl<Integer, ValueType> recDestMap=new ReadOnlyMapImpl<Integer, ValueType>();
private final ReadOnlyMapImpl<Integer, SourceStore> sourceStoreMap=new ReadOnlyMapImpl<Integer, SourceStore>();
private static ParameterMap inputsMap=new ParameterMap();
private static ParameterMap outputsMap=new ParameterMap();
static{
inputsMap.put(StandardParameters.FACTOR);
inputsMap.put(StandardParameters.IN);
outputsMap.put(StandardParameters.OUT);
}
static{
recSourceMap.put(StandardParameters.OUT.i, ValueType.STREAM);
}
private double sampleRate=GlobalParameters.get().getSampleRate();
private class ReconstructionReader implements ProcessElement {
double val;
@Override
public ReadOnlyMap<Integer, ValueType> getSourceTypes() {
return recSourceMap;
}
@Override
public double getValue(int index) {
return val;
}
@Override
public void advanceOutput() {
val=shiftRing[shiftReadPointer%shiftLen];
}
@Override
public void advanceState() {
shiftReadPointer++;
if(shiftReadPointer>shiftLen*1000)
shiftReadPointer-=998*shiftLen;
}
@Override
public String getAbstractName() {
return "";
}
@Override
public String getInstanceName() {
return "";
}
@Override
public void setAbstractName(String abstractName) {
}
@Override
public void setInstanceName(String instanceName) {
}
@Override
public long getNativeSpace() {
return 0;
}
@Override
public ReadOnlyMap<Integer, ValueType> getDestinationTypes() {
return recDestMap;
}
@Override
public void setSource(int connectionIndex, NamedSource source,
int sourceIndex) throws UnknownConnectionException,
TypeMismatchException {
}
@Override
public ReconstructionReader clone() {
return new ReconstructionReader();
}
@Override
public ReadOnlyMap<Integer, SourceStore> getSourceMap() {
return null;
}
/* (non-Javadoc)
* @see de.maramuse.soundcomp.process.ProcessElement#outputsByName()
*/
@Override
public ReadOnlyMap<String, Parameter> outputsByName() {
return null; // not needed, the parser does not see inner elements
}
/* (non-Javadoc)
* @see de.maramuse.soundcomp.process.ProcessElement#inputsByName()
*/
@Override
public ReadOnlyMap<String, Parameter> inputsByName() {
return null; // not needed, the parser does not see inner elements
}
}
private static final ReadOnlyMapImpl<Integer, ValueType> srcMap=new ReadOnlyMapImpl<Integer, ValueType>();
private static final ReadOnlyMapImpl<Integer, ValueType> destMap=new ReadOnlyMapImpl<Integer, ValueType>();
public long nativeSpace=-1;
private ReconstructionReader reconstructionReader=new ReconstructionReader();
static{
srcMap.put(StandardParameters.OUT.i, ValueType.STREAM); // shifted up
destMap.put(StandardParameters.IN.i, ValueType.STREAM);
// the factor by which the frequency is to be multiplied. Must not
// exceed 0,5..2 (an octave)
destMap.put(StandardParameters.FACTOR.i, ValueType.STREAM);
}
StateVariableFilter svf1=new StateVariableFilter();
StateVariableFilter svf2=new StateVariableFilter();
StateVariableFilter reconstruction=new StateVariableFilter();
StateVariableFilter upSampleFilter=new StateVariableFilter();
StateVariableFilter inputHP=new StateVariableFilter();
ConstStream freqHP=ConstStream.c(40); // 40 Hz input HP
ConstStream edgeDetectFilterFreq=ConstStream.c(sampleRate/100);
ConstStream reconstructionFreq=ConstStream.c(sampleRate/20d);
SourceStore factor;
// the analysis ring is a short ring for finding the period in samples
double analysisRing[]=new double[2000];
double inputRing[]=new double[analysisRing.length];
boolean analysisRingFull=false;
int analysisLastIndex=0;
// the shift ring is at upsampled frequency, for resampling.
final int shiftLen=((int)sampleRate);
// the pointer into the shiftring where the current cycle begins
int shiftstart=0;
double shiftRing[]=new double[shiftLen];
int shiftWritePointer=0;
int shiftReadPointer=0;
int analysisPointer=0;
String abstractName, instanceName;
double F, U, U_frac;
int S_F;
double value, input;
final int D=10; // in samples
// the resampleRing should suffice for 0.1 secs even if upsampled by 100, as
// indicated above
double resampleRing[]=new double[100*(int)sampleRate/D];
public TDHS() {
NativeObjects.registerNativeObject(this);
init();
}
TDHS(boolean s) {
init();
}
private void init() {
try{
reconstruction.setSource(StandardParameters.IN.i, reconstructionReader,
StandardParameters.OUT.i);
reconstruction.setSource(StandardParameters.FREQUENCY.i, reconstructionFreq,
StandardParameters.OUT.i);
// no source necessary for reconstructionReader, it is a faked ProcessElement
inputHP.setSource(StandardParameters.FREQUENCY.i, freqHP, StandardParameters.OUT.i);
switch(MODE){
case 1:
svf2.setSource(StandardParameters.IN.i, svf1, StandardParameters.LP.i);
svf2.setSource(StandardParameters.FREQUENCY.i, ConstStream.c(40),
StandardParameters.OUT.i);
svf1.setSource(StandardParameters.FREQUENCY.i, ConstStream.c(40),
StandardParameters.OUT.i);
// the signal for svf1 is set later
break;
case 2:
break;
}
}catch(UnknownConnectionException ex){
// cannot happen, this is constant
}catch(TypeMismatchException ex2){
// cannot happen, this is constant
}
}
@Override
public ReadOnlyMap<Integer, ValueType> getDestinationTypes() {
return destMap;
}
@Override
public void setSource(int connectionIndex, NamedSource source,
int sourceIndex) throws UnknownConnectionException,
TypeMismatchException {
sourceStoreMap.put(connectionIndex, new SourceStore(source, sourceIndex));
if(!source.getSourceTypes().containsKey(sourceIndex))
throw new UnknownConnectionException("attempt to set invalid source for TDHS");
if(connectionIndex==StandardParameters.IN.i){
inputHP.setSource(StandardParameters.IN.i, source, sourceIndex);
if(MODE==1)
svf1.setSource(connectionIndex, source, sourceIndex);
}else if(connectionIndex==StandardParameters.FACTOR.i){
factor=new SourceStore(source, sourceIndex);
}else
throw new UnknownConnectionException("attempt to set unkown TDHS parameter");
}
@Override
public ReadOnlyMap<Integer, ValueType> getSourceTypes() {
return srcMap;
}
@Override
public double getValue(int index) {
return value;
}
@Override
public void advanceOutput() {
reconstruction.advanceOutput();
value=reconstruction.getValue(StandardParameters.LP.i);
}
@Override
public void advanceState() {
boolean periodIsNew=false;
inputHP.advanceState();
inputHP.advanceOutput();
F=factor.getValue();
U=D/F+U_frac;
U_frac=U-Math.floor(U);
U=Math.floor(U);
// if(F<1)System.out.println("pushing "+U+" values with error of "+U_frac);
double _i=inputHP.getValue(StandardParameters.HP.i);
int period=1;
if(MODE!=0){
/*
* here, do some analysis for mode 1 or 2 to determine the period. for certain modes, this MAY depend on a flag
* set by the last determination of the output ring buffer state to avoid unnecessary calculations
*/
switch(MODE){
case 1:
svf1.advanceState();
svf1.advanceOutput();
svf2.advanceState();
svf2.advanceOutput();
double lastVal=analysisRing[(analysisPointer+analysisRing.length-1)%analysisRing.length];
double thisVal=svf2.getValue(StandardParameters.BP.i);
thisVal=_i; // eliminate filters for now
analysisRing[analysisPointer%analysisRing.length]=thisVal;
inputRing[analysisPointer%inputRing.length]=_i;
if(thisVal>0&&lastVal<0){
// rising zero crossing found
period=(analysisPointer-analysisLastIndex);
if(period>analysisRing.length){
period=1; // analysis overflow, no period detectable
}else{
periodIsNew=true;
}
analysisLastIndex=analysisPointer;
}
analysisPointer=analysisPointer+1;
if(analysisPointer>10000*analysisRing.length){
analysisPointer-=9998*analysisRing.length; // avoid numeric wraparaound
}
if(analysisLastIndex>analysisPointer&&analysisLastIndex>9998*analysisRing.length){
analysisLastIndex-=9998*analysisRing.length;
}
break;
case 2:
// get period by covariance
break;
default:
// unknown mode
break;
}
}
/*
* if we're uptuning, we at points have to insert values (the read is too fast). if necessary, do this here by
* duplicating the last period if the modulo difference shiftReadPointer-shiftWritePointer is less than F*two period
*/
if(F>1&&period>1&&periodIsNew){
if(MODE==1){
if(((shiftWritePointer+shiftLen-shiftReadPointer)%shiftLen)<2*D*period){
for(int analysisDuplicator=0; analysisDuplicator<period; analysisDuplicator++){
double v=inputRing[(analysisLastIndex-period+analysisDuplicator)%inputRing.length];
for(int i=1; i<U; i++){
shiftRing[shiftWritePointer++%shiftRing.length]=v;
}
U=D/F+U_frac;
U_frac=U-Math.floor(U);
U=Math.floor(U);
}
}
}else if(MODE==2){
// TODO covariance?
}
}
for(int i=0; i<U; i++){
shiftRing[shiftWritePointer%shiftLen]=_i;
shiftWritePointer++;
}
if(shiftWritePointer>1000*shiftLen){
shiftWritePointer-=998*shiftLen;
shiftReadPointer-=998*shiftLen;
}
/*
* if we're downtuning, we at points have to discard values (write is too fast). if necessary, do this here by
* adjusting the shiftWritePointer backwards by one period if the modulo difference
* shiftReadPointer-shiftWritePointer is larger than F*one period
*/
if(F<1&&periodIsNew&&period>1){
if(MODE==1){
if((shiftWritePointer-shiftReadPointer)>(3*D*period)){
double skip=period*D/F;
int sk=(int)Math.floor(skip);
// System.out.println("have to skip "+period+" input values, going back "+skip+" values");
shiftWritePointer-=sk;
while(shiftWritePointer>shiftReadPointer+(D-1)*period
&&shiftRing[shiftWritePointer%shiftLen]>0)
shiftWritePointer--;
}
}else if(MODE==2){
// TODO covariance?
}
}
for(int i=0; i<D; i++){
reconstructionReader.advanceState();
reconstructionReader.advanceOutput();
reconstruction.advanceState();
}
}
@SuppressWarnings("unused")
private void ignoreForNow() {
// guess a number in the range 80..100 which is a multiple of F
switch((int)Math.floor(10*F)){
case 0: // trim to >= 0.1
S_F=999;
F=0.1;
break;
case 1:
S_F=499;
break;
case 2:
S_F=333;
break;
case 3:
S_F=249;
break;
case 4:
S_F=199;
break;
case 5:
S_F=166;
break; // above should not happen but we will not prevent it
case 6:
S_F=142;
break;
case 7:
S_F=124;
break;
case 8:
S_F=111;
break;
case 9:
S_F=99;
break;
case 10:
S_F=90;
break;
case 11:
S_F=83;
break;
case 12:
S_F=76;
break;
case 13:
S_F=71;
break;
case 14:
S_F=66;
break;
case 15:
S_F=62;
break;
case 16:
S_F=58;
break;
case 17:
S_F=55;
break;
case 18:
S_F=52;
break;
case 19:
S_F=49;
break;
default: // trim to 2
F=2;
S_F=49;
break;
}
input=inputHP.getValue(StandardParameters.HP.i);
shiftRing[shiftWritePointer]=input;
shiftWritePointer=(shiftWritePointer+1)%shiftLen;
analysisRing[analysisPointer]=svf2.getValue(StandardParameters.BP.i);
analysisPointer=(analysisPointer+1)%analysisRing.length;
/*
* As long as we have not filled up the analysisRing at least once, just output zeroes and do nothing but store (we
* do not have the info D to work on yet)
*/
if(!analysisRingFull){ // only the first few milliseconds
if(analysisPointer==0){
analysisRingFull=true;
// D = findPeriod();
shiftstart=0;
shiftWritePointer=-D/2;
}else
return;
}
/*
* factoring scheme: if we shift DOWN (to lower freq), we must REMOVE samples. out of the S samples of a cycle, D
* cycles must get cut off. Therefore, if we reach sample (S-D/2), we start mixing in samples (0-D/2), starting from
* the next Period. This means we always attempt to read at least D/2 behind to be able to mix. When we have reached
* sample D/2, we stop feeding the upsample for another D samples and just store (and get ahead up to 3D/2). Through
* downsampling, the reader will catch up again.
*
* if we shift UP (to higher freq), we must INSERT samples. in addition to a complete cycle, one length D Series is
* inserted (actually from the ringBuf) with overlapping. We also always attempt to read D/2 samples ahead.
*/
if(F>1){
// UPSAMPLING, insert D per S. First find the relevant phase
int relpos=(shiftReadPointer-shiftstart);
if(relpos<0&&relpos>-D/2){
// mix in next period (start-overlap)
}else if(relpos>=0&&relpos<D){
// mix in and out the same period from analysis ring
// this is the actual insertion
}else if(relpos>=D&&relpos<(D*3)/2){
// mix out the rest from analysis ring
}else if(relpos<U-D/2){
// just regular upsample
}else if(relpos<U){
// mix in next period (running-overlap)
}
}else{
// DOWNSAMPLING, remove D out of S
int relpos=(shiftReadPointer-shiftstart);
if(relpos<0&&relpos>-D/2){
// mix in next period (start-overlap)
}else if(relpos==0){
// wait D samples
}else if(relpos<U-D){
// just regular upsample
}else{
// mix in next period (running-overlap)
}
}
/*
* the following only in the OLA windowing regions
*/
// D = findPeriod();
ConstStream upSampled=new ConstStream(input);
try{
upSampleFilter.setSource(StandardParameters.IN.i, upSampled, StandardParameters.OUT.i);
}catch(Exception ex){
// doesn't happen by design
}
}
@Override
public String getAbstractName() {
return abstractName;
}
@Override
public String getInstanceName() {
return instanceName;
}
@Override
public void setAbstractName(String abstractName) {
this.abstractName=abstractName;
}
@Override
public void setInstanceName(String instanceName) {
this.instanceName=instanceName;
}
/**
* @return signal period before upsampling, in samples
*
* the current version just attempts to find consecutive rising zero-crossings. for almost all existing
* signals, this should suffice since the 18dB/oct filter before detection will make the base freq the
* dominating amplitude component.
*
* for future versions, a simplified covariance version is planned. the regular covariance algorithm that is
* supplied with most examples of psola et al. offers much room for optimization, since it contains loop
* invariants in double layered loops.
*
* proposal for an optimized covariance algorithm: - have a "gliding covariance result" for each buffer
* length. it is updated on each cycle by removing the influence of the element leaving the buffer and adding
* the new sample. all these results can be stored in an array of same size as the analysis ring. the
* covariance comparison reduces to one loop: go over all stored previous covariances and adapt them to the
* newest sample as indicated above. while doing so, note the index to the new maximum, so it does not have to
* be searched for. if we store the square of the incoming samples in an additional ring, we need not
* recalculate them for the covariances. so each sample leads to a loop of ring length with only addition,
* subtraction, comparison and array index calculations and therefore is relatively cheap. note that the
* covariance detection is only necessary during the phase starting at (end of cycle minus analysis ring size)
* sample and lasts (analysis ring length) cycles. during the remaining time where we only copy input to
* upsampler, we could skip covariance detection and relevant buffering.
*
*
*/
@SuppressWarnings("unused")
private int findPeriod() {
boolean up=true;
int bl=analysisRing.length;
int firstFound=-1;
double lastVal=analysisRing[(analysisPointer+bl-1)%bl];
if(lastVal>0)
up=false;
for(int i=0; i<analysisRing.length; i++){
// search backwards - newest entries first, to better follow the
// signal
double val=analysisRing[(analysisPointer+bl-i-2)%bl];
if(val*lastVal<0){ // differing sign
if((up&&(val>0))||((!up)&&(val<0))){
if(firstFound==-1)
firstFound=i;
else
return (bl-i+firstFound)%bl;
}
}
}
return bl;
}
@Override
public long getNativeSpace() {
return nativeSpace;
}
@Override
public ReadOnlyMap<Integer, SourceStore> getSourceMap() {
return sourceStoreMap;
}
/**
* @see de.maramuse.soundcomp.process.ProcessElement#clone()
*/
@Override
public TDHS clone() {
TDHS c=new TDHS();
c.abstractName=abstractName;
return c;
}
/*
* (non-Javadoc)
*
* @see de.maramuse.soundcomp.process.ProcessElement#outputsByName()
*/
@Override
public ReadOnlyMap<String, Parameter> outputsByName() {
return outputsMap;
}
/*
* (non-Javadoc)
*
* @see de.maramuse.soundcomp.process.ProcessElement#inputsByName()
*/
@Override
public ReadOnlyMap<String, Parameter> inputsByName() {
return inputsMap;
}
}