Skip to content
Commits on Source (4)
Show whitespace changes
Inline Side-by-side
EEG/EEG.cpp
View file @ 80504b54
......@@ -86,10 +86,9 @@ autoEEG EEG_create (double tmin, double tmax) {
}
integer EEG_getChannelNumber (EEG me, conststring32 channelName) {
for (integer ichan = 1; ichan <= my numberOfChannels; ichan ++) {
for (integer ichan = 1; ichan <= my numberOfChannels; ichan ++)
if (Melder_equ (my channelNames [ichan].get(), channelName))
return ichan;
}
return 0;
}
......@@ -99,7 +98,7 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
char buffer [81];
fread (buffer, 1, 8, f);
buffer [8] = '\0';
bool is24bit = buffer [0] == (char) 255;
const bool is24bit = ( buffer [0] == (char) 255 );
fread (buffer, 1, 80, f);
buffer [80] = '\0';
trace (U"Local subject identification: \"", Melder_peek8to32 (buffer), U"\"");
......@@ -114,27 +113,27 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
trace (U"Start time of recording: \"", Melder_peek8to32 (buffer), U"\"");
fread (buffer, 1, 8, f);
buffer [8] = '\0';
integer numberOfBytesInHeaderRecord = atol (buffer);
const integer numberOfBytesInHeaderRecord = atol (buffer);
trace (U"Number of bytes in header record: ", numberOfBytesInHeaderRecord);
fread (buffer, 1, 44, f);
buffer [44] = '\0';
trace (U"Version of data format: \"", Melder_peek8to32 (buffer), U"\"");
fread (buffer, 1, 8, f);
buffer [8] = '\0';
integer numberOfDataRecords = strtol (buffer, nullptr, 10);
const integer numberOfDataRecords = strtol (buffer, nullptr, 10);
trace (U"Number of data records: ", numberOfDataRecords);
fread (buffer, 1, 8, f);
buffer [8] = '\0';
double durationOfDataRecord = atof (buffer);
const double durationOfDataRecord = atof (buffer);
trace (U"Duration of a data record: ", durationOfDataRecord);
fread (buffer, 1, 4, f);
buffer [4] = '\0';
integer numberOfChannels = atol (buffer);
const integer numberOfChannels = atol (buffer);
trace (U"Number of channels in data record: ", numberOfChannels);
if (numberOfBytesInHeaderRecord != (numberOfChannels + 1) * 256)
Melder_throw (U"Number of bytes in header record (", numberOfBytesInHeaderRecord,
U") doesn't match number of channels (", numberOfChannels, U").");
autostring32vector channelNames (numberOfChannels);
autoSTRVEC channelNames (numberOfChannels);
for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
fread (buffer, 1, 16, f);
buffer [16] = '\0'; // labels of the channels
......@@ -142,16 +141,15 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
* Strip all final spaces.
*/
for (int i = 15; i >= 0; i --) {
if (buffer [i] == ' ') {
if (buffer [i] == ' ')
buffer [i] = '\0';
} else {
else
break;
}
}
channelNames [ichannel] = Melder_8to32 (buffer);
trace (U"Channel <<", channelNames [ichannel].get(), U">>");
}
bool hasLetters = str32equ (channelNames [numberOfChannels].get(), U"EDF Annotations");
const bool hasLetters = str32equ (channelNames [numberOfChannels].get(), U"EDF Annotations");
double samplingFrequency = undefined;
for (integer channel = 1; channel <= numberOfChannels; channel ++) {
fread (buffer, 1, 80, f);
......@@ -193,7 +191,7 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
for (integer channel = 1; channel <= numberOfChannels; channel ++) {
fread (buffer, 1, 8, f);
buffer [8] = '\0'; // number of samples in each data record
integer numberOfSamplesInThisDataRecord = atol (buffer);
const integer numberOfSamplesInThisDataRecord = atol (buffer);
if (isundef (samplingFrequency)) {
numberOfSamplesPerDataRecord = numberOfSamplesInThisDataRecord;
samplingFrequency = numberOfSamplesInThisDataRecord / durationOfDataRecord;
......@@ -207,36 +205,37 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
fread (buffer, 1, 32, f);
buffer [32] = '\0'; // reserved
}
double duration = numberOfDataRecords * durationOfDataRecord;
const double duration = numberOfDataRecords * durationOfDataRecord;
autoEEG him = EEG_create (0, duration);
his numberOfChannels = numberOfChannels;
autoSound me = Sound_createSimple (numberOfChannels, duration, samplingFrequency);
Melder_assert (my nx == numberOfSamplesPerDataRecord * numberOfDataRecords);
autoNUMvector <uint8> dataBuffer ((integer) 0, 3 * numberOfSamplesPerDataRecord - 1);
autoBYTEVEC dataBuffer = newBYTEVECzero (3 * numberOfSamplesPerDataRecord);
for (integer record = 1; record <= numberOfDataRecords; record ++) {
for (integer channel = 1; channel <= numberOfChannels; channel ++) {
double factor = channel == numberOfChannels ? 1.0 : physicalMinimum [channel] / digitalMinimum [channel];
if (channel < numberOfChannels - EEG_getNumberOfExtraSensors (him.get())) factor /= 1000000.0;
double factor = ( channel == numberOfChannels ? 1.0 : physicalMinimum [channel] / digitalMinimum [channel] );
if (channel < numberOfChannels - EEG_getNumberOfExtraSensors (him.get()))
factor /= 1000000.0;
if (is24bit) {
fread (& dataBuffer [0], 3, (size_t) numberOfSamplesPerDataRecord, f);
uint8 *p = & dataBuffer [0];
fread (dataBuffer.asArgumentToFunctionThatExpectsZeroBasedArray(), 3, (size_t) numberOfSamplesPerDataRecord, f);
byte *p = & dataBuffer [1];
for (integer i = 1; i <= numberOfSamplesPerDataRecord; i ++) {
integer sample = i + (record - 1) * numberOfSamplesPerDataRecord;
const integer sample = i + (record - 1) * numberOfSamplesPerDataRecord;
Melder_assert (sample <= my nx);
uint8 lowByte = *p ++, midByte = *p ++, highByte = *p ++;
const uint8 lowByte = *p ++, midByte = *p ++, highByte = *p ++;
uint32 externalValue = ((uint32) highByte << 16) | ((uint32) midByte << 8) | (uint32) lowByte;
if ((highByte & 128) != 0) // is the 24-bit sign bit on?
externalValue |= 0xFF00'0000; // extend negative sign to 32 bits
my z [channel] [sample] = (int32) externalValue * factor;
}
} else {
fread (& dataBuffer [0], 2, (size_t) numberOfSamplesPerDataRecord, f);
uint8 *p = & dataBuffer [0];
fread (dataBuffer.asArgumentToFunctionThatExpectsZeroBasedArray(), 2, (size_t) numberOfSamplesPerDataRecord, f);
byte *p = & dataBuffer [1];
for (integer i = 1; i <= numberOfSamplesPerDataRecord; i ++) {
integer sample = i + (record - 1) * numberOfSamplesPerDataRecord;
const integer sample = i + (record - 1) * numberOfSamplesPerDataRecord;
Melder_assert (sample <= my nx);
uint8 lowByte = *p ++, highByte = *p ++;
uint16 externalValue = (uint16) ((uint16) highByte << 8) | (uint16) lowByte;
const uint8 lowByte = *p ++, highByte = *p ++;
const uint16 externalValue = (uint16) ((uint16) highByte << 8) | (uint16) lowByte;
my z [channel] [sample] = (int16) externalValue * factor;
}
}
......@@ -244,21 +243,20 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
}
int numberOfStatusBits = 8;
for (integer i = 1; i <= my nx; i ++) {
uint32 value = (uint32) (int32) my z [numberOfChannels] [i];
if (value & 0x0000'FF00) {
const uint32 value = (uint32) (int32) my z [numberOfChannels] [i];
if (value & 0x0000'FF00)
numberOfStatusBits = 16;
}
}
autoTextGrid thee;
if (hasLetters) {
thee = TextGrid_create (0, duration, U"Mark Trigger", U"Mark Trigger");
autoMelderString letters;
double time = undefined;
for (integer i = 1; i <= my nx; i ++) {
uint32 value = (uint32) (int32) my z [numberOfChannels] [i];
const uint32 value = (uint32) (int32) my z [numberOfChannels] [i];
for (int ibyte = 1; ibyte <= numberOfStatusBits / 8; ibyte ++) {
uint32 mask = ( ibyte == 1 ? 0x0000'00ff : 0x0000'ff00 );
char32 kar = ( ibyte == 1 ? (value & mask) : (value & mask) >> 8 );
const uint32 mask = ( ibyte == 1 ? 0x0000'00ff : 0x0000'ff00 );
const char32 kar = ( ibyte == 1 ? (value & mask) : (value & mask) >> 8 );
if (kar != U'\0' && kar != 20) {
MelderString_appendCharacter (& letters, kar);
} else if (letters. string [0] != U'\0') {
......@@ -307,13 +305,13 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
U""
);
for (int bit = 1; bit <= numberOfStatusBits; bit ++) {
uint32 bitValue = 1 << (bit - 1);
const uint32 bitValue = 1 << (bit - 1);
IntervalTier tier = (IntervalTier) thy tiers->at [bit];
for (integer i = 1; i <= my nx; i ++) {
uint32 previousValue = i == 1 ? 0 : (uint32) (int32) my z [numberOfChannels] [i - 1];
uint32 thisValue = (uint32) (int32) my z [numberOfChannels] [i];
const uint32 previousValue = ( i == 1 ? 0 : (uint32) (int32) my z [numberOfChannels] [i - 1] );
const uint32 thisValue = (uint32) (int32) my z [numberOfChannels] [i];
if ((thisValue & bitValue) != (previousValue & bitValue)) {
double time = i == 1 ? 0.0 : my x1 + (i - 1.5) * my dx;
const double time = ( i == 1 ? 0.0 : my x1 + (i - 1.5) * my dx );
if (time != 0.0)
TextGrid_insertBoundary (thee.get(), bit, time);
if ((thisValue & bitValue) != 0)
......@@ -432,7 +430,7 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
}
static void detrend (VEC const& channel) {
double firstValue = channel [1], lastValue = channel [channel.size];
const double firstValue = channel [1], lastValue = channel [channel.size];
channel [1] = channel [channel.size] = 0.0;
for (integer isamp = 2; isamp < channel.size; isamp ++)
channel [isamp] -= ((isamp - 1.0) * lastValue + (channel.size - isamp) * firstValue) / (channel.size - 1);
......@@ -529,7 +527,7 @@ void EEG_setChannelToZero (EEG me, integer channelNumber) {
void EEG_setChannelToZero (EEG me, conststring32 channelName) {
try {
integer channelNumber = EEG_getChannelNumber (me, channelName);
const integer channelNumber = EEG_getChannelNumber (me, channelName);
if (channelNumber == 0)
Melder_throw (U"No channel named \"", channelName, U"\".");
EEG_setChannelToZero (me, channelNumber);
......@@ -554,7 +552,7 @@ autoEEG EEG_extractChannel (EEG me, integer channelNumber) {
Melder_throw (U"No channel ", channelNumber, U".");
autoEEG thee = EEG_create (my xmin, my xmax);
thy numberOfChannels = 1;
thy channelNames = autostring32vector (1);
thy channelNames = autoSTRVEC (1);
thy channelNames [1] = Melder_dup (my channelNames [1].get());
thy sound = Sound_extractChannel (my sound.get(), channelNumber);
thy textgrid = Data_copy (my textgrid.get());
......@@ -566,7 +564,7 @@ autoEEG EEG_extractChannel (EEG me, integer channelNumber) {
autoEEG EEG_extractChannel (EEG me, conststring32 channelName) {
try {
integer channelNumber = EEG_getChannelNumber (me, channelName);
const integer channelNumber = EEG_getChannelNumber (me, channelName);
if (channelNumber == 0)
Melder_throw (U"No channel named \"", channelName, U"\".");
return EEG_extractChannel (me, channelNumber);
......@@ -577,15 +575,15 @@ autoEEG EEG_extractChannel (EEG me, conststring32 channelName) {
autoEEG EEG_extractChannels (EEG me, constVECVU const& channelNumbers) {
try {
integer numberOfChannels = channelNumbers.size;
const integer numberOfChannels = channelNumbers.size;
Melder_require (numberOfChannels > 0,
U"The number of channels should be greater than 0.");
autoEEG you = EEG_create (my xmin, my xmax);
your sound = Sound_extractChannels (my sound.get(), channelNumbers);
your numberOfChannels = numberOfChannels;
your channelNames = autostring32vector (numberOfChannels);
your channelNames = autoSTRVEC (numberOfChannels);
for (integer ichan = 1; ichan <= numberOfChannels; ichan ++) {
integer originalChannelNumber = Melder_iround (channelNumbers [ichan]);
const integer originalChannelNumber = Melder_iround (channelNumbers [ichan]);
your channelNames [ichan] = Melder_dup (my channelNames [originalChannelNumber].get());
}
your textgrid = Data_copy (my textgrid.get());
......@@ -614,9 +612,8 @@ void EEG_removeChannel (EEG me, integer channelNumber) {
try {
if (channelNumber < 1 || channelNumber > my numberOfChannels)
Melder_throw (U"No channel ", channelNumber, U".");
for (integer ichan = channelNumber; ichan < my numberOfChannels; ichan ++) {
for (integer ichan = channelNumber; ichan < my numberOfChannels; ichan ++)
my channelNames [ichan] = my channelNames [ichan + 1].move();
}
my channelNames [my numberOfChannels]. reset();
my numberOfChannels -= 1;
Sound_removeChannel (my sound.get(), channelNumber);
......@@ -627,7 +624,7 @@ void EEG_removeChannel (EEG me, integer channelNumber) {
void EEG_removeChannel (EEG me, conststring32 channelName) {
try {
integer channelNumber = EEG_getChannelNumber (me, channelName);
const integer channelNumber = EEG_getChannelNumber (me, channelName);
if (channelNumber == 0)
Melder_throw (U"No channel named \"", channelName, U"\".");
EEG_removeChannel (me, channelNumber);
......@@ -641,8 +638,8 @@ autoEEG EEGs_concatenate (OrderedOf<structEEG>* me) {
if (my size < 1)
Melder_throw (U"Cannot concatenate zero EEG objects.");
EEG first = my at [1];
integer numberOfChannels = first -> numberOfChannels;
autostring32vector channelNames = newSTRVECcopy (first -> channelNames.get());
const integer numberOfChannels = first -> numberOfChannels;
autoSTRVEC channelNames = newSTRVECcopy (first -> channelNames.get());
for (integer ieeg = 2; ieeg <= my size; ieeg ++) {
EEG other = my at [ieeg];
if (other -> numberOfChannels != numberOfChannels)
......
FFNet/FFNet.cpp
View file @ 80504b54
......@@ -77,7 +77,8 @@ autostring32 FFNet_createNameFromTopology (FFNet me) {
static double sigmoid (FFNet /*me*/, double x, double *out_deriv) {
const double act = NUMsigmoid (x);
if (out_deriv) *out_deriv = act * (1.0 - act);
if (out_deriv)
*out_deriv = act * (1.0 - act);
return act;
}
......@@ -111,7 +112,6 @@ static double minimumCrossEntropy (FFNet me, constVEC& target) {
for (integer i = 1; i <= my numberOfOutputs; i ++, k ++) {
const double t1 = 1.0 - target [i];
const double o1 = 1.0 - my activity [k];
cost -= target [i] * log (my activity [k]) + t1 * log (o1);
my error [k] = -t1 / o1 + target [i] / my activity [k];
}
......@@ -134,8 +134,9 @@ static void bookkeeping (FFNet me) {
U"Number of weights is incorrect.");
my numberOfWeights = numberOfWeights;
// The following test is essential because when an FFNet is read from file the w array already exists
/*
The following test is essential because when an FFNet is read from file the w array already exists
*/
if (! my w.at)
my w = newVECzero (my numberOfWeights);
......@@ -257,7 +258,8 @@ void FFNet_setCostFunction (FFNet me, int costType) {
double FFNet_getBias (FFNet me, integer layer, integer unit) {
try {
const integer node = FFNet_getNodeNumberFromUnitNumber (me, unit, layer);
Melder_require (node > 0, U"Not a valid unit / layer combination.");
Melder_require (node > 0,
U"Not a valid unit / layer combination.");
const integer bias_unit = my wLast [node];
return my w [bias_unit];
} catch (MelderError) {
......@@ -312,7 +314,7 @@ void FFNet_reset (FFNet me, double weightRange) {
conststring32 FFNet_getCategoryOfOutputUnit (FFNet me, integer outputUnit) {
conststring32 result = U"-- undefined --";
if (my outputCategories && outputUnit <= my outputCategories -> size) {
SimpleString ss = my outputCategories->at [outputUnit];
const SimpleString ss = my outputCategories->at [outputUnit];
result = ss -> string.get();
}
return result;
......@@ -322,7 +324,7 @@ integer FFNet_getOutputUnitOfCategory (FFNet me, const char32* category) {
integer result = 0;
if (my outputCategories) {
for (integer i = 1; i <= my outputCategories -> size; i ++) {
SimpleString s = my outputCategories->at [i];
const SimpleString s = my outputCategories->at [i];
if (Melder_equ (s -> string.get(), category)) {
result = i;
break;
......@@ -336,11 +338,15 @@ integer FFNet_getOutputUnitOfCategory (FFNet me, const char32* category) {
/* step 1 */
void FFNet_propagate (FFNet me, constVEC input, autoVEC *output) {
Melder_assert (my numberOfInputs == input.size);
// clamp input pattern on the network
/*
Clamp input pattern on the network
*/
my activity.part (1, my numberOfInputs) <<= input;
// on hidden units use activation function
integer k = 1, numberOfNodes = my outputsAreLinear ? my numberOfNodes - my numberOfOutputs : my numberOfNodes;
/*
On hidden units use activation function
*/
const integer numberOfNodes = my outputsAreLinear ? my numberOfNodes - my numberOfOutputs : my numberOfNodes;
integer k = 1;
for (integer i = my numberOfInputs + 2; i <= numberOfNodes; i ++) {
if (my isbias [i])
continue;
......@@ -350,7 +356,9 @@ void FFNet_propagate (FFNet me, constVEC input, autoVEC *output) {
my activity [i] = my nonLinearity (me, act, & my deriv [i]);
}
// on output units use another activation function
/*
On output units use another activation function
*/
if (my outputsAreLinear) {
for (integer i = numberOfNodes + 1; i <= my numberOfNodes; i ++) {
if (my isbias [i])
......@@ -373,12 +381,15 @@ void FFNet_propagate (FFNet me, constVEC input, autoVEC *output) {
double FFNet_computeError (FFNet me, constVEC target) {
Melder_assert (my numberOfOutputs == target.size);
// compute error at output layer
/*
Compute error at output layer
*/
const double cost = my costFunction (me, target);
for (integer i = 1; i <= my numberOfNodes - my numberOfOutputs; i ++)
my error [i] = 0.0;
// backpropagation of errors from output to first hidden layer
/*
Backpropagation of errors from output to first hidden layer
*/
for (integer i = my numberOfNodes; i > my numberOfInputs + 1; i--) {
if (my isbias [i])
continue;
......@@ -395,12 +406,11 @@ double FFNet_computeError (FFNet me, constVEC target) {
void FFNet_computeDerivative (FFNet me) {
integer k = 1;
for (integer i = my numberOfInputs + 2; i <= my numberOfNodes; i ++) {
for (integer i = my numberOfInputs + 2; i <= my numberOfNodes; i ++)
if (! my isbias [i])
for (integer node = my nodeFirst [i]; node <= my nodeLast [i]; node ++, k ++)
my dwi [k] = - my error [i] * my activity [node];
}
}
/******* end operation ******************************************************/
......@@ -413,10 +423,9 @@ integer FFNet_getWinningUnit (FFNet me, integer labeling) {
sum += my activity [k + ioutput];
const double random = NUMrandomUniform (0.0, sum);
for (winningUnit = my numberOfOutputs; winningUnit >= 2; winningUnit--) {
for (winningUnit = my numberOfOutputs; winningUnit >= 2; winningUnit--)
if (random > (sum -= my activity [k + winningUnit]))
break;
}
} else { /* winner-takes-all */
double max = my activity [k + 1];
for (integer ioutput = 2; ioutput <= my numberOfOutputs; ioutput ++)
......@@ -546,7 +555,9 @@ void FFNet_drawTopology (FFNet me, Graphics g) {
const double y2WC = dy / 2 + layer * dy;
double dx2 = dx, x2WC;
double x2 = (maxNumOfUnits - numberOfUnitsInLayer + 1) * dx2 / 2;
/* draw the units */
/*
Draw the units
*/
if (! dxIsFixed) {
dx2 = 1.0 / numberOfUnitsInLayer;
x2 = dx2 / 2.0;
......@@ -629,8 +640,8 @@ void FFNet_drawActivation (FFNet me, Graphics g) {
}
x2WC = x2;
for (integer iunit = 1; iunit <= numberOfUnitsInLayer; iunit ++, node ++) {
double activity = my activity [node];
double radius = r1 * (fabs (activity) < 0.05 ? 0.05 : fabs (activity));
const double activity = my activity [node];
const double radius = r1 * (fabs (activity) < 0.05 ? 0.05 : fabs (activity));
/*Graphics_setColour (g, activity < 0 ? Melder_BLACK : Melder_RED);*/
Graphics_circle (g, x2WC, y2WC, radius);
if (activity < 0)
......
FFNet/FFNet_ActivationList_Categories.cpp
View file @ 80504b54
......@@ -28,10 +28,9 @@
static integer winnerTakesAll (FFNet me, constVEC activation) {
integer pos = 1;
double max = activation[1];
for (integer i = 2; i <= my numberOfOutputs; i ++) {
for (integer i = 2; i <= my numberOfOutputs; i ++)
if (activation [i] > max)
max = activation [pos = i];
}
return pos;
}
......@@ -43,7 +42,8 @@ static integer stochastic (FFNet me, constVEC activation) {
const double number = NUMrandomUniform (0.0, range);
for (i = 1; i <= my numberOfOutputs; i ++) {
lower += activation [i];
if (number < lower) break;
if (number < lower)
break;
}
return i;
}
......@@ -74,17 +74,16 @@ autoActivationList FFNet_Categories_to_ActivationList (FFNet me, Categories thee
autoCategories uniq = Categories_selectUniqueItems (thee);
Melder_require (my outputCategories,
U"The FFNet does not have categories.");
const integer nl = OrderedOfString_isSubsetOf (uniq.get(), my outputCategories.get(), 0);
Melder_require (nl > 0,
U"The Categories should match the categories of the FFNet.");
autoActivationList him = ActivationList_create (thy size, my numberOfOutputs);
for (integer i = 1; i <= thy size; i ++) {
SimpleString category = thy at [i];
const SimpleString category = thy at [i];
const integer pos = OrderedOfString_indexOfItem_c (my outputCategories.get(), category -> string.get());
if (pos < 1)
Melder_throw (U"The FFNet doesn't know the category ", category -> string.get(), U".");
Melder_require (pos > 0,
U"The FFNet doesn't know the category ", category -> string.get(), U".");
his z [i] [pos] = 1.0;
}
return him;
......
FFNet/FFNet_Eigen.cpp
View file @ 80504b54
......@@ -50,28 +50,26 @@ void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pc
Graphics_setWindow (g, xmin, xmax, ymin, ymax);
for (integer i = 1; i <= my numberOfUnitsInLayer [1]; i ++) {
const integer unitOffset = my numberOfInputs + 1;
const double bias = my w [my wLast [unitOffset + i]];
double c1 = 0.0, c2 = 0.0;
double x [6], y [6], xs [3], ys [3];
integer ns = 0;
for (integer j = 1; j <= my numberOfInputs; j ++) {
c1 += my w [my wFirst [unitOffset + i] + j - 1] * eigen -> eigenvectors [ix] [j];
c2 += my w [my wFirst [unitOffset + i] + j - 1] * eigen -> eigenvectors [iy] [j];
}
double x [6], y [6], xs [3], ys [3];
x [1] = x [2] = x [5] = xmin;
x [3] = x [4] = xmax;
y [1] = y [4] = y [5] = ymin;
y [2] = y [3] = ymax;
integer ns = 0;
const double bias = my w [my wLast [unitOffset + i]];
for (integer j = 1; j <= 4; j++) {
const double p1 = c1 * x [j] + c2 * y [j] + bias;
const double p2 = c1 * x [j + 1] + c2 * y [j + 1] + bias;
const double r = fabs (p1) / (fabs (p1) + fabs (p2));
if (p1 *p2 > 0 || r == 0.0)
continue;
if (++ ns > 2)
break;
xs [ns] = x [j] + (x [j + 1] - x [j]) * r;
ys [ns] = y [j] + (y [j + 1] - y [j]) * r;
}
......@@ -88,8 +86,7 @@ void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pc
from layer j with the plane spanned by eigenvectors pcx and pcy.
*/
void FFNet_Eigen_drawDecisionPlaneInEigenspace (FFNet me, Eigen thee, Graphics g, integer unit, integer layer,
integer pcx, integer pcy, double xmin, double xmax, double ymin, double ymax)
{
integer pcx, integer pcy, double xmin, double xmax, double ymin, double ymax) {
if (layer < 1 || layer > my numberOfLayers)
return;
if (unit < 1 || unit > my numberOfUnitsInLayer [layer])
......@@ -100,7 +97,6 @@ void FFNet_Eigen_drawDecisionPlaneInEigenspace (FFNet me, Eigen thee, Graphics g
if (numberOfUnitsInLayer_m1 != thy dimension)
return;
double x1, x2, y1, y2;
Graphics_inqWindow (g, & x1, & x2, & y1, & y2);
if (xmax <= xmin) {
......@@ -169,7 +165,7 @@ void FFNet_Eigen_drawDecisionPlaneInEigenspace (FFNet me, Eigen thee, Graphics g
"for unit ", unit, U" in layer ", layer, U" with the plane spanned by the eigenvectors because \nboth planes are parallel.");
return;
}
double xi [3], yi [3]; /* Intersections */
double xi [3], yi [3]; // Intersections
const double ni = NUMgetIntersectionsWithRectangle (x1, y1, x2, y2, xmin, ymin, xmax, ymax, xi, yi);
if (ni == 2)
Graphics_line (g, xi [1], yi [1], xi [2], yi [2]);
......
FFNet/FFNet_PatternList_ActivationList.cpp
View file @ 80504b54
......@@ -29,21 +29,22 @@
static double func (Daata object, VEC const& p) {
FFNet me = (FFNet) object;
Minimizer thee = my minimizer.get();
longdouble fp = 0.0;
const Minimizer thee = my minimizer.get();
for (integer j = 1, k = 1; k <= my numberOfWeights; k ++) {
my dw [k] = 0.0;
if (my wSelected [k])
my w [k] = p [j ++];
}
longdouble fp = 0.0;
for (integer i = 1; i <= my numberOfPatterns; i ++) {
FFNet_propagate (me, my inputPattern.row (i), nullptr);
fp += FFNet_computeError (me, my targetActivation.row (i));
FFNet_computeDerivative (me);
/* derivative (cumulative) */
for (integer k = 1; k <= my numberOfWeights; k ++)
my dw [k] += my dwi [k];
/*
Derivative (cumulative)
*/
my dw.part (1, my numberOfWeights) += my dwi.part (1, my numberOfWeights);
}
thy numberOfFunctionCalls ++;
return (double) fp;
......@@ -66,19 +67,18 @@ static void _FFNet_PatternList_ActivationList_checkDimensions (FFNet me, Pattern
U"The Activation and the FFNet do not match.\nThe number of columns in the Activation must equal the number of outputs in the FFNet.");
Melder_require (p -> ny == a -> ny,
U"The PatternList and the ActivationList do not match.\nThe number of rows in the PatternList must equal the number of rows in the Activation.");
if (! _PatternList_checkElements (p)) {
Melder_throw (U"All PatternList elements should be in the interval [0, 1].\nYou could use \"Formula...\" to scale the PatternList values first."); }
if (! _ActivationList_checkElements (a)) {
Melder_throw (U"All Activation elements should be in the interval [0, 1].\nYou could use \"Formula...\" to scale the Activation values first.");
}
Melder_require (_PatternList_checkElements (p),
U"All PatternList elements should be in the interval [0, 1].\nYou could use \"Formula...\" to scale the PatternList values first.");
Melder_require (_ActivationList_checkElements (a),
U"All Activation elements should be in the interval [0, 1].\nYou could use \"Formula...\" to scale the Activation values first.");
}
static void _FFNet_PatternList_ActivationList_learn (FFNet me, PatternList pattern, ActivationList activation, integer maxNumOfEpochs, double tolerance, int costFunctionType, bool reset) {
try {
_FFNet_PatternList_ActivationList_checkDimensions (me, pattern, activation);
// Link the things to be learned
/*
Link the things to be learned
*/
my numberOfPatterns = pattern -> ny;
my inputPattern = pattern -> z.get();
my targetActivation = activation -> z.get();
......@@ -87,11 +87,9 @@ static void _FFNet_PatternList_ActivationList_learn (FFNet me, PatternList patte
if (reset) {
autoVEC wbuf = newVECzero (my dimension);
integer k = 1;
for (integer i = 1; i <= my numberOfWeights; i ++) {
if (my wSelected [i]) {
for (integer i = 1; i <= my numberOfWeights; i ++)
if (my wSelected [i])
wbuf [k ++] = my w [i];
}
}
Minimizer_reset (my minimizer.get(), wbuf.get());
}
......@@ -112,12 +110,16 @@ static void _FFNet_PatternList_ActivationList_learn (FFNet me, PatternList patte
void FFNet_PatternList_ActivationList_learnSD (FFNet me, PatternList p, ActivationList a, integer maxNumOfEpochs, double tolerance, double learningRate, double momentum, int costFunctionType) {
bool resetMinimizer = false;
/* Did we choose another minimizer */
/*
Did we choose another minimizer
*/
if (my minimizer && ! Thing_isa (my minimizer.get(), classSteepestDescentMinimizer)) {
my minimizer.reset();
resetMinimizer = true;
}
/* create the minimizer if it doesn't exist */
/*
Create the minimizer if it doesn't exist
*/
if (! my minimizer) {
resetMinimizer = true;
my minimizer = SteepestDescentMinimizer_create (my dimension, me, func, dfunc_optimized);
......@@ -129,14 +131,16 @@ void FFNet_PatternList_ActivationList_learnSD (FFNet me, PatternList p, Activati
void FFNet_PatternList_ActivationList_learnSM (FFNet me, PatternList p, ActivationList a, integer maxNumOfEpochs, double tolerance, int costFunctionType) {
bool resetMinimizer = false;
// Did we choose another minimizer
/*
Did we choose another minimizer
*/
if (my minimizer.get() && ! Thing_isa (my minimizer.get(), classVDSmagtMinimizer)) {
my minimizer.reset();
resetMinimizer = true;
}
// create the minimizer if it doesn't exist
/*
Create the minimizer if it doesn't exist
*/
if (! my minimizer.get()) {
resetMinimizer = true;
my minimizer = VDSmagtMinimizer_create (my dimension, me, func, dfunc_optimized);
......@@ -176,10 +180,8 @@ autoActivationList FFNet_PatternList_to_ActivationList (FFNet me, PatternList p,
const integer numberOfPatterns = p -> ny;
autoActivationList thee = ActivationList_create (numberOfPatterns, my numberOfUnitsInLayer [layer]);
for (integer i = 1; i <= numberOfPatterns; i ++) {
for (integer i = 1; i <= numberOfPatterns; i ++)
FFNet_propagateToLayer (me, p -> z.row (i), thy z.row (i), layer);
}
return thee;
} catch (MelderError) {
Melder_throw (me, U": no ActivationList created.");
......
LPC/Cepstrogram.cpp
View file @ 80504b54
......@@ -66,37 +66,33 @@ autoPowerCepstrogram PowerCepstrogram_create (double tmin, double tmax, integer
void PowerCepstrogram_paint (PowerCepstrogram me, Graphics g, double tmin, double tmax, double qmin, double qmax, double dBmaximum, int autoscaling, double dynamicRangedB, double dynamicCompression, bool garnish) {
Function_unidirectionalAutowindow (me, & tmin, & tmax);
if (qmax <= qmin) { qmin = my ymin; qmax = my ymax; }
if (qmax <= qmin) {
qmin = my ymin;
qmax = my ymax;
}
integer itmin, itmax, ifmin, ifmax;
if (! Matrix_getWindowSamplesX (me, tmin - 0.49999 * my dx, tmax + 0.49999 * my dx, & itmin, & itmax) ||
! Matrix_getWindowSamplesY (me, qmin - 0.49999 * my dy, qmax + 0.49999 * my dy, & ifmin, & ifmax)) {
if (Matrix_getWindowSamplesX (me, tmin - 0.49999 * my dx, tmax + 0.49999 * my dx, & itmin, & itmax) == 0 ||
Matrix_getWindowSamplesY (me, qmin - 0.49999 * my dy, qmax + 0.49999 * my dy, & ifmin, & ifmax) == 0)
return;
}
autoMatrix thee = Data_copy (me);
double min = 1e308, max = -min;
for (integer i = 1; i <= my ny; i ++) {
for (integer j = 1; j <= my nx; j ++) {
double val = TO10LOG (my z [i] [j]);
if (val < min) min = val;
if (val > max) max = val;
thy z [i] [j] = val;
MelderExtremaWithInit extrema;
for (integer irow = 1; irow <= my ny; irow ++) {
for (integer icol = 1; icol <= my nx; icol ++) {
double val = TO10LOG (my z [irow] [icol]);
extrema.update (val);
thy z [irow] [icol] = val;
}
}
double dBminimum = dBmaximum - dynamicRangedB;
if (autoscaling) {
dBminimum = min;
dBmaximum = max;
dBminimum = extrema.min;
dBmaximum = extrema.max;
}
for (integer j = 1; j <= my nx; j ++) {
double lmax = thy z [1] [j];
for (integer i = 2; i <= my ny; i ++) {
if (thy z [i] [j] > lmax) {
lmax = thy z [i] [j];
}
}
double factor = dynamicCompression * (max - lmax);
thy z.column(j) += factor;
for (integer icol = 1; icol <= my nx; icol ++) {
const double lmax = NUMmax (thy z.column (icol));
const double factor = dynamicCompression * (extrema.max - lmax);
thy z.column(icol) += factor;
}
Graphics_setInner (g);
......@@ -121,10 +117,10 @@ void PowerCepstrogram_paint (PowerCepstrogram me, Graphics g, double tmin, doubl
void PowerCepstrogram_subtractTrend_inplace (PowerCepstrogram me, double qstartFit, double qendFit, kCepstrumTrendType lineType, kCepstrumTrendFit fitMethod) {
try {
autoPowerCepstrum thee = PowerCepstrum_create (my ymax, my ny);
for (integer i = 1; i <= my nx; i ++) {
thy z.row (1) <<= my z.column (i);
for (integer icol = 1; icol <= my nx; icol ++) {
thy z.row (1) <<= my z.column (icol);
PowerCepstrum_subtractTrend_inplace (thee.get(), qstartFit, qendFit, lineType, fitMethod);
my z.column (i) <<= thy z.row (1);
my z.column (icol) <<= thy z.row (1);
}
} catch (MelderError) {
Melder_throw (me, U": no tilt subtracted (inline).");
......@@ -146,14 +142,15 @@ autoTable PowerCepstrogram_to_Table_hillenbrand (PowerCepstrogram me, double pit
try {
autoTable thee = Table_createWithColumnNames (my nx, U"time quefrency cpp f0");
autoPowerCepstrum him = PowerCepstrum_create (my ymax, my ny);
for (integer i = 1; i <= my nx; i ++) {
his z.row (1) <<= my z.column (i);
double qpeak, cpp = PowerCepstrum_getPeakProminence_hillenbrand (him.get(), pitchFloor, pitchCeiling, & qpeak);
double time = Sampled_indexToX (me, i);
Table_setNumericValue (thee.get(), i, 1, time);
Table_setNumericValue (thee.get(), i, 2, qpeak);
Table_setNumericValue (thee.get(), i, 3, cpp); // Cepstrogram_getCPPS depends on this index 3!!
Table_setNumericValue (thee.get(), i, 4, 1.0 / qpeak);
for (integer icol = 1; icol <= my nx; icol ++) {
his z.row (1) <<= my z.column (icol);
double qpeak;
const double cpp = PowerCepstrum_getPeakProminence_hillenbrand (him.get(), pitchFloor, pitchCeiling, & qpeak);
const double time = Sampled_indexToX (me, icol);
Table_setNumericValue (thee.get(), icol, 1, time);
Table_setNumericValue (thee.get(), icol, 2, qpeak);
Table_setNumericValue (thee.get(), icol, 3, cpp); // Cepstrogram_getCPPS depends on this index 3!!
Table_setNumericValue (thee.get(), icol, 4, 1.0 / qpeak);
}
return thee;
} catch (MelderError) {
......@@ -165,17 +162,17 @@ autoTable PowerCepstrogram_to_Table_cpp (PowerCepstrogram me, double pitchFloor,
try {
autoTable thee = Table_createWithColumnNames (my nx, U"time quefrency cpp f0 rnr");
autoPowerCepstrum him = PowerCepstrum_create (my ymax, my ny);
for (integer i = 1; i <= my nx; i ++) {
his z.row (1) <<= my z.column (i);
for (integer icol = 1; icol <= my nx; icol ++) {
his z.row (1) <<= my z.column (icol);
double qpeak, cpp = PowerCepstrum_getPeakProminence (him.get(), pitchFloor, pitchCeiling, interpolation,
qstartFit, qendFit, lineType, fitMethod, & qpeak);
double rnr = PowerCepstrum_getRNR (him.get(), pitchFloor, pitchCeiling, deltaF0);
double time = Sampled_indexToX (me, i);
Table_setNumericValue (thee.get(), i, 1, time);
Table_setNumericValue (thee.get(), i, 2, qpeak);
Table_setNumericValue (thee.get(), i, 3, cpp); // Cepstrogram_getCPPS depends on this index!!
Table_setNumericValue (thee.get(), i, 4, 1.0 / qpeak);
Table_setNumericValue (thee.get(), i, 5, rnr);
double time = Sampled_indexToX (me, icol);
Table_setNumericValue (thee.get(), icol, 1, time);
Table_setNumericValue (thee.get(), icol, 2, qpeak);
Table_setNumericValue (thee.get(), icol, 3, cpp); // Cepstrogram_getCPPS depends on this index!!
Table_setNumericValue (thee.get(), icol, 4, 1.0 / qpeak);
Table_setNumericValue (thee.get(), icol, 5, rnr);
}
return thee;
} catch (MelderError) {
......@@ -186,7 +183,9 @@ autoTable PowerCepstrogram_to_Table_cpp (PowerCepstrogram me, double pitchFloor,
autoPowerCepstrogram PowerCepstrogram_smooth (PowerCepstrogram me, double timeAveragingWindow, double quefrencyAveragingWindow) {
try {
autoPowerCepstrogram thee = Data_copy (me);
// 1. average across time
/*
1. average across time
*/
integer numberOfFrames = Melder_ifloor (timeAveragingWindow / my dx);
if (numberOfFrames > 1) {
autoVEC qin = newVECraw (my nx);
......@@ -195,7 +194,9 @@ autoPowerCepstrogram PowerCepstrogram_smooth (PowerCepstrogram me, double timeAv
VECsmoothByMovingAverage_preallocated (thy z.row (iq), qin.get(), numberOfFrames);
}
}
// 2. average across quefrencies
/*
2. average across quefrencies
*/
integer numberOfQuefrencyBins = Melder_ifloor (quefrencyAveragingWindow / my dy);
if (numberOfQuefrencyBins > 1) {
autoVEC qin = newVECraw (thy ny);
......@@ -244,32 +245,34 @@ autoPowerCepstrogram Matrix_to_PowerCepstrogram (Matrix me) {
autoPowerCepstrogram Sound_to_PowerCepstrogram (Sound me, double pitchFloor, double dt, double maximumFrequency, double preEmphasisFrequency) {
try {
// minimum analysis window has 3 periods of lowest pitch
double analysisWidth = 3.0 / pitchFloor;
double windowDuration = 2.0 * analysisWidth; /* gaussian window */
double analysisWidth = 3.0 / pitchFloor; // minimum analysis window has 3 periods of lowest pitch
double windowDuration = 2.0 * analysisWidth; // gaussian window
integer nFrames;
// Convenience: analyse the whole sound into one Cepstrogram_frame
if (windowDuration > my dx * my nx) {
if (windowDuration > my dx * my nx)
windowDuration = my dx * my nx;
}
double t1, samplingFrequency = 2 * maximumFrequency;
const double samplingFrequency = 2 * maximumFrequency;
autoSound sound = Sound_resample (me, samplingFrequency, 50);
Sound_preEmphasis (sound.get(), preEmphasisFrequency);
double t1;
Sampled_shortTermAnalysis (me, windowDuration, dt, & nFrames, & t1);
autoSound sframe = Sound_createSimple (1, windowDuration, samplingFrequency);
autoSound window = Sound_createGaussian (windowDuration, samplingFrequency);
// find out the size of the FFT
/*
Find out the size of the FFT
*/
integer nfft = 2;
while (nfft < sframe -> nx) nfft *= 2;
integer nq = nfft / 2 + 1;
double qmax = 0.5 * nfft / samplingFrequency, dq = qmax / (nq - 1);
while (nfft < sframe -> nx)
nfft *= 2;
const integer nq = nfft / 2 + 1;
const double qmax = 0.5 * nfft / samplingFrequency, dq = qmax / (nq - 1);
autoPowerCepstrogram thee = PowerCepstrogram_create (my xmin, my xmax, nFrames, dt, t1, 0, qmax, nq, dq, 0);
autoMelderProgress progress (U"Cepstrogram analysis");
for (integer iframe = 1; iframe <= nFrames; iframe++) {
double t = Sampled_indexToX (thee.get(), iframe);
const double t = Sampled_indexToX (thee.get(), iframe);
Sound_into_Sound (sound.get(), sframe.get(), t - windowDuration / 2);
Vector_subtractMean (sframe.get());
Sounds_multiply (sframe.get(), window.get());
......@@ -295,7 +298,7 @@ autoCepstrum Spectrum_to_Cepstrum_hillenbrand (Spectrum me) {
// originalNumberOfSamplesProbablyOdd irrelevant
Melder_require (my x1 == 0.0,
U"A Fourier-transformable Spectrum should have a first frequency of 0 Hz, not ", my x1, U" Hz.");
integer numberOfSamples = my nx - 1;
const integer numberOfSamples = my nx - 1;
autoCepstrum thee = Cepstrum_create (0.5 / my dx, my nx);
NUMfft_Table_init (& fftTable, my nx);
autoVEC amp = newVECraw (my nx);
......@@ -321,10 +324,10 @@ autoCepstrum Spectrum_to_Cepstrum_hillenbrand (Spectrum me) {
static void complexfftoutput_to_power (constVEC fft, VEC dbs, bool to_db) {
double valsq = fft [1] * fft [1];
dbs[1] = to_db ? TOLOG (valsq) : valsq;
integer nfftdiv2p1 = (fft.size + 2) / 2;
integer nend = fft.size % 2 == 0 ? nfftdiv2p1 : nfftdiv2p1 + 1;
const integer nfftdiv2p1 = (fft.size + 2) / 2;
const integer nend = fft.size % 2 == 0 ? nfftdiv2p1 : nfftdiv2p1 + 1;
for (integer i = 2; i < nend; i ++) {
double re = fft [i + i - 2], im = fft [i + i - 1];
const double re = fft [i + i - 2], im = fft [i + i - 1];
valsq = re * re + im * im;
dbs [i] = to_db ? TOLOG (valsq) : valsq;
}
......@@ -348,22 +351,26 @@ autoPowerCepstrogram Sound_to_PowerCepstrogram_hillenbrand (Sound me, double pit
} else {
thee = Data_copy (me);
}
// pre-emphasis with fixed coefficient 0.9
/*
Pre-emphasis with fixed coefficient 0.9
*/
for (integer i = thy nx; i > 1; i --)
thy z [1] [i] -= 0.9 * thy z [1] [i - 1];
integer nosInWindow = Melder_ifloor (analysisWidth * samplingFrequency), numberOfFrames;
const integer nosInWindow = Melder_ifloor (analysisWidth * samplingFrequency);
Melder_require (nosInWindow >= 8,
U"Analysis window too short.");
double t1;
integer numberOfFrames;
Sampled_shortTermAnalysis (thee.get(), analysisWidth, dt, & numberOfFrames, & t1);
autoVEC hamming = newVECraw (nosInWindow);
for (integer i = 1; i <= nosInWindow; i ++)
hamming [i] = 0.54 - 0.46 * cos (NUM2pi * (i - 1) / (nosInWindow - 1));
integer nfft = 8; // minimum possible
while (nfft < nosInWindow) { nfft *= 2; }
while (nfft < nosInWindow)
nfft *= 2;
const integer nfftdiv2 = nfft / 2;
autoVEC fftbuf = newVECzero (nfft); // "complex" array
autoVEC spectrum = newVECzero (nfftdiv2 + 1); // +1 needed
......@@ -377,8 +384,10 @@ autoPowerCepstrogram Sound_to_PowerCepstrogram_hillenbrand (Sound me, double pit
for (integer iframe = 1; iframe <= numberOfFrames; iframe ++) {
const double tbegin = std::max (thy xmin, t1 + (iframe - 1) * dt - analysisWidth / 2.0);
const integer istart = std::max (integer (1), Sampled_xToLowIndex (thee.get(), tbegin)); // ppgb: afronding naar beneden?
const integer iend = istart + nosInWindow - 1; // ppgb: crash if not provably iend <= thy nx
const integer istart = std::max (1_integer, Sampled_xToLowIndex (thee.get(), tbegin)); // ppgb: afronding naar beneden?
integer iend = istart + nosInWindow - 1;
if (iend > thy nx)
iend = thy nx;
fftbuf.part (1, nosInWindow) <<= thy z.row (1).part (istart, iend) * hamming.all();
fftbuf.part (nosInWindow + 1, nfft) <<= 0.0;
......@@ -421,7 +430,7 @@ double PowerCepstrogram_getCPPS (PowerCepstrogram me, bool subtractTiltBeforeSmo
autoPowerCepstrogram smooth = PowerCepstrogram_smooth (flattened ? flattened.get() : me, timeAveragingWindow, quefrencyAveragingWindow);
autoTable table = PowerCepstrogram_to_Table_cpp (smooth.get(), pitchFloor, pitchCeiling, deltaF0, interpolation, qstartFit, qendFit, lineType, fitMethod);
double cpps = Table_getMean (table.get(), 3);
const double cpps = Table_getMean (table.get(), 3);
return cpps;
} catch (MelderError) {
Melder_throw (me, U": no CPPS value calculated.");
......@@ -436,7 +445,7 @@ double PowerCepstrogram_getCPPS_hillenbrand (PowerCepstrogram me, bool subtractT
autoPowerCepstrogram smooth = PowerCepstrogram_smooth (subtractTiltBeforeSmoothing ? him.get() : me, timeAveragingWindow, quefrencyAveragingWindow);
autoTable table = PowerCepstrogram_to_Table_hillenbrand (smooth.get(), pitchFloor, pitchCeiling);
double cpps = Table_getMean (table.get(), 3);
const double cpps = Table_getMean (table.get(), 3);
return cpps;
} catch (MelderError) {
Melder_throw (me, U": no CPPS value calculated.");
......
LPC/Cepstrum.cpp
View file @ 80504b54
......@@ -97,15 +97,13 @@ static void _Cepstrum_draw (Cepstrum me, Graphics g, double qmin, double qmax, d
}
integer imin, imax;
if (! Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax)) {
if (Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax) == 0)
return;
}
integer numberOfSelected = imax - imin + 1;
autoVEC y = newVECraw (numberOfSelected);
for (integer i = 1; i <= numberOfSelected; i ++) {
for (integer i = 1; i <= numberOfSelected; i ++)
y [i] = my v_getValueAtSample (imin + i - 1, (power ? 1 : 0), 0);
}
if (autoscaling)
NUMextrema (y.get(), & minimum, & maximum);
......@@ -145,9 +143,8 @@ void PowerCepstrum_drawTrendLine (PowerCepstrum me, Graphics g, double qmin, dou
if (dBminimum >= dBmaximum) { // autoscaling
integer imin, imax;
if (! Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax)) {
if (Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax) == 0)
return;
}
integer numberOfPoints = imax - imin + 1;
dBminimum = dBmaximum = my v_getValueAtSample (imin, 1, 0);
for (integer i = 2; i <= numberOfPoints; i ++) {
......@@ -219,7 +216,7 @@ void PowerCepstrum_fitTrendLine (PowerCepstrum me, double qmin, double qmax, dou
}
integer imin, imax;
if (! Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax))
if (Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax) == 0)
return;
if (imin == 1 && lineType == kCepstrumTrendType::ExponentialDecay)
imin = 2; // because log(0) is undefined
......@@ -345,7 +342,7 @@ double PowerCepstrum_getRNR (PowerCepstrum me, double pitchFloor, double pitchCe
double qmin = 1.0 / pitchCeiling, qmax = 1.0 / pitchFloor, peakdB, qpeak;
PowerCepstrum_getMaximumAndQuefrency (me, pitchFloor, pitchCeiling, 2, & peakdB, & qpeak);
integer imin, imax;
if (! Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax))
if (Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax) == 0)
return rnr;
integer ndata = imax - imin + 1;
......
LPC/Cepstrum_and_Spectrum.cpp
View file @ 80504b54
......@@ -32,7 +32,7 @@
autoPowerCepstrum Spectrum_to_PowerCepstrum (Spectrum me) {
try {
autoSpectrum dBspectrum = Data_copy (me);
VEC re = dBspectrum -> z.row (1), im = dBspectrum -> z.row (2);
const VEC re = dBspectrum -> z.row (1), im = dBspectrum -> z.row (2);
for (integer i = 1; i <= dBspectrum -> nx; i ++) {
re [i] = log (re [i] * re [i] + im [i] * im [i] + 1e-300);
im [i] = 0.0;
......@@ -40,7 +40,7 @@ autoPowerCepstrum Spectrum_to_PowerCepstrum (Spectrum me) {
autoSound cepstrum = Spectrum_to_Sound (dBspectrum.get());
autoPowerCepstrum thee = PowerCepstrum_create (0.5 / my dx, my nx);
for (integer i = 1; i <= thy nx; i ++) {
double val = cepstrum -> z [1] [i];
const double val = cepstrum -> z [1] [i];
thy z [1] [i] = val * val;
}
return thee;
......@@ -52,17 +52,14 @@ autoPowerCepstrum Spectrum_to_PowerCepstrum (Spectrum me) {
autoCepstrum Spectrum_to_Cepstrum (Spectrum me) {
try {
autoSpectrum dBspectrum = Data_copy (me);
VEC re = dBspectrum -> z.row (1), im = dBspectrum -> z.row (2);
const VEC re = dBspectrum -> z.row (1), im = dBspectrum -> z.row (2);
for (integer i = 1; i <= dBspectrum -> nx; i ++) {
re [i] = log (re [i] * re [i] + im [i] * im [i] + 1e-300);
im [i] = 0.0;
}
autoSound cepstrum = Spectrum_to_Sound (dBspectrum.get());
autoCepstrum thee = Cepstrum_create (0.5 / my dx, my nx);
for (integer i = 1; i <= thy nx; i ++) {
double val = cepstrum -> z [1] [i];
thy z [1] [i] = val;
}
thy z.row (1) <<= cepstrum -> z.row(1).part (1, thy nx);
return thee;
} catch (MelderError) {
Melder_throw (me, U": not converted to Sound.");
......@@ -73,12 +70,11 @@ autoSpectrum Cepstrum_to_Spectrum (Cepstrum me) { //TODO power cepstrum
try {
autoCepstrum cepstrum = Data_copy (me);
cepstrum -> z [1] [1] = my z [1] [1];
for (integer i = 2; i <= cepstrum -> nx; i ++) {
cepstrum -> z [1] [i] = 2 * my z [1] [i];
}
for (integer i = 2; i <= cepstrum -> nx; i ++)
cepstrum -> z [1] [i] = 2.0 * my z [1] [i];
autoSpectrum thee = Sound_to_Spectrum ((Sound) cepstrum.get(), true);
VEC re = thy z.row (1), im = thy z.row (2);
const VEC re = thy z.row (1), im = thy z.row (2);
for (integer i = 1; i <= thy nx; i ++) {
re [i] = exp (0.5 * re [i]); // i.e., sqrt (exp(re [i]))
im [i] = 0.0;
......
LPC/Cepstrumc.cpp
View file @ 80504b54
......@@ -63,16 +63,14 @@ void Cepstrumc_Frame_init (Cepstrumc_Frame me, int nCoefficients) {
my nCoefficients = nCoefficients;
}
void Cepstrumc_init (Cepstrumc me, double tmin, double tmax, integer nt, double dt, double t1,
int nCoefficients, double samplingFrequency) {
void Cepstrumc_init (Cepstrumc me, double tmin, double tmax, integer nt, double dt, double t1, integer nCoefficients, double samplingFrequency) {
my samplingFrequency = samplingFrequency;
my maxnCoefficients = nCoefficients;
Sampled_init (me, tmin, tmax, nt, dt, t1);
my frame = newvectorzero <structCepstrumc_Frame> (nt);
}
autoCepstrumc Cepstrumc_create (double tmin, double tmax, integer nt, double dt, double t1,
int nCoefficients, double samplingFrequency) {
autoCepstrumc Cepstrumc_create (double tmin, double tmax, integer nt, double dt, double t1, integer nCoefficients, double samplingFrequency) {
try {
autoCepstrumc me = Thing_new (Cepstrumc);
Cepstrumc_init (me.get(), tmin, tmax, nt, dt, t1, nCoefficients, samplingFrequency);
......@@ -82,88 +80,91 @@ autoCepstrumc Cepstrumc_create (double tmin, double tmax, integer nt, double dt,
}
}
static void regression (Cepstrumc me, integer frame, double r [], integer nr) {
integer nc = 1e6; double sumsq = 0;
for (integer i = 0; i <= my maxnCoefficients; i ++) {
r [i] = 0;
}
if (frame <= nr / 2 || frame >= my nx - nr / 2) {
static void regression (VEC const& r, Cepstrumc me, integer frameNumber, integer numberOfRegressionFrames) {
Melder_assert (r.size == my maxnCoefficients + 1);
r <<= 0.0;
if (frameNumber <= numberOfRegressionFrames / 2 || frameNumber > my nx - numberOfRegressionFrames / 2)
return;
}
for (integer j = -nr / 2; j <= nr / 2; j ++) {
Cepstrumc_Frame f = & my frame [frame + j];
if (f->nCoefficients < nc) {
integer nc = INTEGER_MAX;
longdouble sumsq = 0.0;
for (integer j = -numberOfRegressionFrames / 2; j <= numberOfRegressionFrames / 2; j ++) {
const Cepstrumc_Frame f = & my frame [frameNumber + j];
if (f -> nCoefficients < nc)
nc = f -> nCoefficients;
}
sumsq += j * j;
}
for (integer i = 0; i <= nc; i ++) {
for (integer j = -nr / 2; j <= nr / 2; j ++) {
Cepstrumc_Frame f = & my frame [frame + j];
r [i] += f->c [i] * j / sumsq / my dx;
for (integer j = -numberOfRegressionFrames / 2; j <= numberOfRegressionFrames / 2; j ++) {
const Cepstrumc_Frame f = & my frame [frameNumber + j];
r [i + 1] += f -> c [i] * j / (double (sumsq) * my dx);
}
}
}
autoDTW Cepstrumc_to_DTW (Cepstrumc me, Cepstrumc thee, double wc, double wle, double wr, double wer, double dtr, int matchStart, int matchEnd, int constraint) {
try {
integer nr = Melder_ifloor (dtr / my dx);
if (my maxnCoefficients != thy maxnCoefficients) {
Melder_throw (U"Cepstrumc orders should be equal.");
}
if (wr != 0 && nr < 2) {
Melder_throw (U"Time window for regression coefficients too small.");
}
if (nr % 2 == 0) {
nr ++;
}
if (wr != 0) {
Melder_casual (U"Number of frames used for regression coefficients ", nr);
}
integer numberOfRegressionFrames = Melder_ifloor (dtr / my dx);
Melder_require (my maxnCoefficients == thy maxnCoefficients,
U"Cepstrumc orders should be equal.");
Melder_require (! (wr != 0.0 && numberOfRegressionFrames < 2),
U"Time window for regression coefficients too small.");
if (numberOfRegressionFrames % 2 == 0)
numberOfRegressionFrames += 1;
if (wr != 0.0)
Melder_casual (U"Number of frames used for regression coefficients ", numberOfRegressionFrames);
autoDTW him = DTW_create (my xmin, my xmax, my nx, my dx, my x1, thy xmin, thy xmax, thy nx, thy dx, thy x1);
autoNUMvector <double> ri ((integer) 0, my maxnCoefficients);
autoNUMvector <double> rj ((integer) 0, my maxnCoefficients);
// Calculate distance matrix
autoVEC ri = newVECraw (my maxnCoefficients + 1);
autoVEC rj = newVECraw (my maxnCoefficients + 1);
/*
Calculate distance matrix.
*/
autoMelderProgress progress (U"");
for (integer i = 1; i <= my nx; i ++) {
Cepstrumc_Frame fi = & my frame [i];
regression (me, i, ri.peek(), nr);
for (integer j = 1; j <= thy nx; j ++) {
Cepstrumc_Frame fj = & thy frame [j];
longdouble d, dist = 0.0, distr = 0.0;
if (wc != 0) { /* cepstral distance */
for (integer iframe = 1; iframe <= my nx; iframe ++) {
const Cepstrumc_Frame fi = & my frame [iframe];
regression (ri.get(), me, iframe, numberOfRegressionFrames);
for (integer jframe = 1; jframe <= thy nx; jframe ++) {
const Cepstrumc_Frame fj = & thy frame [jframe];
longdouble dist = 0.0, distr = 0.0;
/*
Cepstral distance.
*/
if (wc != 0) {
for (integer k = 1; k <= fj -> nCoefficients; k ++) {
d = fi -> c [k] - fj -> c [k];
dist += d * d;
const double dci = fi -> c [k] - fj -> c [k];
dist += dci * dci;
}
dist *= wc;
}
// log energy distance
d = fi -> c [0] - fj -> c [0];
dist += wle * d * d;
if (wr != 0) { // regression distance
regression (thee, j, rj.peek(), nr);
/*
Log energy distance.
*/
const double dc0 = fi -> c [0] - fj -> c [0];
dist += wle * dc0 * dc0;
/*
Regression distance.
*/
if (wr != 0) {
regression (rj.get(), thee, jframe, numberOfRegressionFrames);
for (integer k = 1; k <= fj -> nCoefficients; k ++) {
d = ri [k] - rj [k];
distr += d * d;
const double drci = ri [k + 1] - rj [k + 1];
distr += drci * drci;
}
dist += wr * distr;
}
if (wer != 0) { // regression on c [0]: log(energy)
if (wr == 0) {
regression (thee, j, rj.peek(), nr);
}
d = ri [0] - rj [0];
dist += wer * d * d;
/*
Regression on c [0]: log(energy)
*/
if (wer != 0) {
if (wr == 0)
regression (rj.get(), thee, jframe, numberOfRegressionFrames);
const double drc0 = ri [1] - rj [1];
dist += wer * drc0 * drc0;
}
dist /= wc + wle + wr + wer;
his z [i] [j] = sqrt ((double) dist); // prototype along y-direction
his z [iframe] [jframe] = sqrt ((double) dist); // prototype along y-direction
}
Melder_progress ( (double) i / my nx, U"Calculate distances: frame ",
i, U" from ", my nx, U".");
Melder_progress ( (double) iframe / my nx, U"Calculate distances: frame ", iframe, U" from ", my nx, U".");
}
DTW_findPath (him.get(), matchStart, matchEnd, constraint);
return him;
......@@ -174,15 +175,13 @@ autoDTW Cepstrumc_to_DTW (Cepstrumc me, Cepstrumc thee, double wc, double wle, d
autoMatrix Cepstrumc_to_Matrix (Cepstrumc me) {
try {
autoMatrix thee = Matrix_create (my xmin, my xmax, my nx, my dx, my x1,
0, my maxnCoefficients, my maxnCoefficients + 1, 1, 0);
autoMatrix thee = Matrix_create (my xmin, my xmax, my nx, my dx, my x1, 0, my maxnCoefficients, my maxnCoefficients + 1, 1, 0);
for (integer i = 1; i <= my nx; i ++) {
Cepstrumc_Frame him = & my frame [i];
for (integer j = 1; j <= his nCoefficients + 1; j ++) {
const Cepstrumc_Frame him = & my frame [i];
for (integer j = 1; j <= his nCoefficients + 1; j ++)
thy z [j] [i] = his c [j - 1];
}
}
return thee;
} catch (MelderError) {
Melder_throw (me, U": no Matrix created.");
......
LPC/Cepstrumc.h
View file @ 80504b54
......@@ -32,10 +32,10 @@
#include "Cepstrumc_def.h"
void Cepstrumc_init (Cepstrumc me, double tmin, double tmax, integer nt, double dt, double t1,
int nCoefficients, double samplingFrequency);
integer nCoefficients, double samplingFrequency);
autoCepstrumc Cepstrumc_create (double tmin, double tmax, integer nt, double dt, double t1,
int nCoefficients, double samplingFrequency);
integer nCoefficients, double samplingFrequency);
/******************* Frames ************************************************/
......
LPC/Formant_extensions.cpp
View file @ 80504b54
/* Formant_extensions.cpp
*
* Copyright (C) 2012-2017 David Weenink
* Copyright (C) 2012-2019 David Weenink
*
* This code is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -22,37 +22,42 @@
void Formant_formula (Formant me, double tmin, double tmax, integer formantmin, integer formantmax, Interpreter interpreter, conststring32 expression) {
try {
Function_unidirectionalAutowindow (me, & tmin, & tmax);
integer numberOfPossibleFormants = my maxnFormants;
const integer numberOfPossibleFormants = my maxnFormants;
if (formantmax >= formantmin) {
formantmin = 1;
formantmax = numberOfPossibleFormants;
}
Melder_clipLeft (integer (1), & formantmin);
Melder_clipLeft (1_integer, & formantmin);
Melder_clipRight (& formantmax, numberOfPossibleFormants);
autoMatrix fb = Matrix_create (my xmin, my xmax, my nx, my dx, my x1, 1.0, 2 * numberOfPossibleFormants, 2 * numberOfPossibleFormants, 1.0, 1.0);
for (integer iframe = 1; iframe <= my nx; iframe ++) {
Formant_Frame frame = & my frames [iframe];
integer numberOfFormants = std::min (integer (frame -> nFormants), numberOfPossibleFormants);
for (integer iformant = 1; iformant <= numberOfFormants; iformant++) {
const Formant_Frame frame = & my frames [iframe];
const integer numberOfFormants = std::min (integer (frame -> nFormants), numberOfPossibleFormants);
for (integer iformant = 1; iformant <= numberOfFormants; iformant++)
if (iformant <= frame -> nFormants) {
fb -> z [2 * iformant - 1] [iframe] = frame -> formant [iformant]. frequency;
fb -> z [2 * iformant ] [iframe] = frame -> formant [iformant]. bandwidth;
}
}
}
// Apply formula
/*
Apply the formaula
*/
const double ymin = 2.0 * formantmin - 1.0, ymax = 2.0 * formantmax;
Matrix_formula_part (fb.get(), tmin, tmax, ymin, ymax, expression, interpreter, nullptr);
// Put results back in Formant
/*
Put results back in Formant
*/
integer ixmin, ixmax, iymin, iymax;
(void) Matrix_getWindowSamplesX (fb.get(), tmin, tmax, & ixmin, & ixmax);
(void) Matrix_getWindowSamplesY (fb.get(), ymin, ymax, & iymin, & iymax);
for (integer iframe = ixmin; iframe <= ixmax; iframe ++) {
// if some of the formant frequencies are set to zero => remove the formant
Formant_Frame frame = & my frames [iframe];
integer numberOfFormants = std::min (integer (frame -> nFormants), formantmax);
/*
If some of the formant frequencies are set to zero => remove the formant
*/
const Formant_Frame frame = & my frames [iframe];
const integer numberOfFormants = std::min (integer (frame -> nFormants), formantmax);
integer iformantto = ( formantmin > 1 ? formantmin - 1 : 0 );
for (integer iformant = formantmin; iformant <= numberOfFormants; iformant++) {
const double frequency = fb -> z [2 * iformant - 1] [iframe];
......@@ -61,22 +66,22 @@ void Formant_formula (Formant me, double tmin, double tmax, integer formantmin,
iformantto ++;
frame -> formant [iformantto]. frequency = frequency;
frame -> formant [iformantto]. bandwidth = bandWidth;
} else {
} else
frame -> formant [iformant]. frequency = frame -> formant [iformant]. bandwidth = 0.0;
}
}
// shift the (higher) formants down if necessary.
/*
Shift the (higher) formants down if necessary.
*/
for (integer iformant = formantmax + 1; iformant <= frame -> nFormants; iformant ++) {
double frequency = fb -> z [2 * iformant - 1] [iframe];
double bandWidth = fb -> z [2 * iformant ] [iframe];
const double frequency = fb -> z [2 * iformant - 1] [iframe];
const double bandWidth = fb -> z [2 * iformant ] [iframe];
if (frequency > 0 && bandWidth > 0) {
iformantto ++;
frame -> formant [iformantto]. frequency = frequency;
frame -> formant [iformantto]. bandwidth = bandWidth;
} else {
} else
frame -> formant [iformant]. frequency = frame -> formant [iformant]. bandwidth = 0.0;
}
}
frame -> nFormants = iformantto;
}
} catch (MelderError) {
......@@ -86,15 +91,16 @@ void Formant_formula (Formant me, double tmin, double tmax, integer formantmin,
autoIntensityTier Formant_Spectrogram_to_IntensityTier (Formant me, Spectrogram thee, integer iformant) {
try {
Melder_require (my xmin == thy xmin && my xmax == thy xmax, U"The start and end times of the Formant and the Spectrogram should be equal.");
Melder_require (iformant > 0 && iformant <= my maxnFormants, U"Formant number should be in the range [1, ", my maxnFormants, U"].");
Melder_require (my xmin == thy xmin && my xmax == thy xmax,
U"The start and end times of the Formant and the Spectrogram should be equal.");
Melder_require (iformant > 0 && iformant <= my maxnFormants,
U"Formant number should be in the range [1, ", my maxnFormants, U"].");
autoIntensityTier him = IntensityTier_create (my xmin, my xmax);
double previousValue = -80000.0; // can never occur
double previousTime = my xmin;
for (integer iframe = 1; iframe <= my nx; iframe ++) {
Formant_Frame frame = & my frames [iframe];
integer numberOfFormants = frame -> nFormants;
const Formant_Frame frame = & my frames [iframe];
const integer numberOfFormants = frame -> nFormants;
const double time = Sampled_indexToX (me, iframe);
double value = 0.0;
if (iformant <= numberOfFormants) {
......
LPC/LPC.cpp
View file @ 80504b54
/* LPC.cpp
*
* Copyright (C) 1994-2018 David Weenink
* Copyright (C) 1994-2019 David Weenink
*
* This code is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -89,7 +89,7 @@ void LPC_drawGain (LPC me, Graphics g, double tmin, double tmax, double gmin, do
if (! Sampled_getWindowSamples (me, tmin, tmax, & itmin, & itmax))
return;
integer numberOfSelected = itmax - itmin + 1;
const integer numberOfSelected = itmax - itmin + 1;
autoVEC gain = newVECraw (numberOfSelected);
for (integer iframe = itmin; iframe <= itmax; iframe ++)
......@@ -106,7 +106,7 @@ void LPC_drawGain (LPC me, Graphics g, double tmin, double tmax, double gmin, do
Graphics_setInner (g);
Graphics_setWindow (g, tmin, tmax, gmin, gmax);
for (integer iframe = itmin; iframe <= itmax; iframe ++) {
double x = Sampled_indexToX (me, iframe);
const double x = Sampled_indexToX (me, iframe);
Graphics_speckle (g, x, gain [iframe - itmin + 1]);
}
Graphics_unsetInner (g);
......@@ -129,7 +129,7 @@ autoMatrix LPC_downto_Matrix_lpc (LPC me) {
try {
autoMatrix thee = Matrix_create (my xmin, my xmax, my nx, my dx, my x1, 0.5, 0.5 + my maxnCoefficients, my maxnCoefficients, 1.0, 1.0);
for (integer j = 1; j <= my nx; j ++) {
LPC_Frame lpc = & my d_frames [j];
const LPC_Frame lpc = & my d_frames [j];
thy z.column (j) <<= lpc-> a.get();
}
return thee;
......@@ -143,7 +143,7 @@ autoMatrix LPC_downto_Matrix_rc (LPC me) {
autoMatrix thee = Matrix_create (my xmin, my xmax, my nx, my dx, my x1, 0.5, 0.5 + my maxnCoefficients, my maxnCoefficients, 1.0, 1.0);
autoVEC rc = newVECzero (my maxnCoefficients);
for (integer j = 1; j <= my nx; j ++) {
LPC_Frame lpc = & my d_frames [j];
const LPC_Frame lpc = & my d_frames [j];
VECrc_from_lpc (rc.part (1, lpc -> nCoefficients), lpc -> a.part (1, lpc -> nCoefficients));
if (lpc -> nCoefficients < my maxnCoefficients)
rc.part (lpc -> nCoefficients + 1, my maxnCoefficients) <<= 0.0;
......@@ -161,7 +161,7 @@ autoMatrix LPC_downto_Matrix_area (LPC me) {
autoVEC rc = newVECraw (my maxnCoefficients);
autoVEC area = newVECraw (my maxnCoefficients);
for (integer j = 1; j <= my nx; j ++) {
LPC_Frame lpc = & my d_frames [j];
const LPC_Frame lpc = & my d_frames [j];
VECrc_from_lpc (rc.part (1, lpc -> nCoefficients), lpc -> a.part (1, lpc -> nCoefficients));
VECarea_from_rc (area.part (1, lpc -> nCoefficients), rc.part (1, lpc -> nCoefficients));
if (lpc -> nCoefficients < my maxnCoefficients)
......
LPC/LPC_and_Cepstrumc.cpp
View file @ 80504b54
/* LPC_and_Cepstrumc.cpp
*
* Copyright (C) 1994-2018 David Weenink
* Copyright (C) 1994-2019 David Weenink
*
* This code is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -24,14 +24,15 @@
#include "LPC_and_Cepstrumc.h"
void LPC_Frame_into_Cepstrumc_Frame (LPC_Frame me, Cepstrumc_Frame thee) {
integer n = std::min (my nCoefficients, thy nCoefficients);
double *c = thy c, *a = my a.at;
const integer n = std::min (my nCoefficients, thy nCoefficients);
const double *a = my a.at;
double *c = thy c;
c [0] = 0.5 * log (my gain);
if (n == 0)
return;
c [1] = -a [1];
for (integer i = 2; i <= n; i ++) {
c [i] = 0;
c [i] = 0.0;
for (integer k = 1; k < i; k ++)
c [i] += a [i - k] * c [k] * k;
c [i] = -a [i] - c [i] / i;
......@@ -52,6 +53,9 @@ void Cepstrumc_Frame_into_LPC_Frame (Cepstrumc_Frame me, LPC_Frame thee) {
a [i] += a [j] * c [i - j];
a [i] /= -i;
}
/*
Undo the modification of the c array
*/
for (integer i = 2; i <= thy nCoefficients; i ++)
c [i] /= i;
}
......@@ -71,8 +75,7 @@ autoCepstrumc LPC_to_Cepstrumc (LPC me) {
autoLPC Cepstrumc_to_LPC (Cepstrumc me) {
try {
autoLPC thee = LPC_create (my xmin, my xmax, my nx, my dx, my x1,
my maxnCoefficients, 1.0 / my samplingFrequency);
autoLPC thee = LPC_create (my xmin, my xmax, my nx, my dx, my x1, my maxnCoefficients, 1.0 / my samplingFrequency);
for (integer i = 1; i <= my nx; i ++) {
LPC_Frame_init (& thy d_frames [i], my frame [i].nCoefficients);
Cepstrumc_Frame_into_LPC_Frame (& my frame [i], & thy d_frames [i]);
......
LPC/LPC_and_Formant.cpp
View file @ 80504b54
/* LPC_and_Formant.cpp
*
* Copyright (C) 1994-2018 David Weenink
* Copyright (C) 1994-2019 David Weenink
*
* This code is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -28,10 +28,9 @@
void Formant_Frame_init (Formant_Frame me, integer numberOfFormants) {
my nFormants = numberOfFormants;
if (numberOfFormants > 0) {
if (numberOfFormants > 0)
my formant = newvectorzero <structFormant_Formant> (my nFormants);
}
}
void Formant_Frame_scale (Formant_Frame me, double scale) {
for (integer i = 1; i <= my nFormants; i ++) {
......@@ -41,21 +40,19 @@ void Formant_Frame_scale (Formant_Frame me, double scale) {
}
void Roots_into_Formant_Frame (Roots me, Formant_Frame thee, double samplingFrequency, double margin) {
integer n = my max - my min + 1;
const integer n = my max - my min + 1;
autoVEC fc = newVECzero (n);
autoVEC bc = newVECzero (n);
// Determine the formants and bandwidths
/*
Determine the formants and bandwidths
*/
thy nFormants = 0;
double fLow = margin, fHigh = samplingFrequency / 2 - margin;
const double fLow = margin, fHigh = samplingFrequency / 2 - margin;
for (integer i = my min; i <= my max; i ++) {
if (my v [i].im < 0) {
if (my v [i].im < 0)
continue;
}
double f = fabs (atan2 (my v [i].im, my v [i].re)) * samplingFrequency / 2.0 / NUMpi;
const double f = fabs (atan2 (my v [i].im, my v [i].re)) * samplingFrequency / 2.0 / NUMpi;
if (f >= fLow && f <= fHigh) {
/*b = - log (my v [i].re * my v [i].re + my v [i].im * my v [i].im) * samplingFrequency / 2 / NUMpi;*/
double b = - log (dcomplex_abs (my v [i])) * samplingFrequency / NUMpi;
thy nFormants ++;
fc [thy nFormants] = f;
......@@ -84,9 +81,9 @@ void LPC_Frame_into_Formant_Frame (LPC_Frame me, Formant_Frame thee, double samp
autoFormant LPC_to_Formant (LPC me, double margin) {
try {
const double samplingFrequency = 1.0 / my samplingPeriod;
integer nmax = my maxnCoefficients;
const integer nmax = my maxnCoefficients;
integer numberOfSuspectFrames = 0;
integer interval = ( nmax > 20 ? 1 : 10 );
const integer interval = ( nmax > 20 ? 1 : 10 );
Melder_require (nmax < 100,
U"We cannot find the roots of a polynomial of order > 99.");
Melder_require (margin < samplingFrequency / 4.0,
......@@ -96,20 +93,20 @@ autoFormant LPC_to_Formant (LPC me, double margin) {
autoMelderProgress progress (U"LPC to Formant");
for (integer i = 1; i <= my nx; i ++) {
Formant_Frame formant = & thy frames [i];
LPC_Frame lpc = & my d_frames [i];
// Initialisation of Formant_Frame is taken care of in Roots_into_Formant_Frame!
for (integer iframe = 1; iframe <= my nx; iframe ++) {
const Formant_Frame formant = & thy frames [iframe];
const LPC_Frame lpc = & my d_frames [iframe];
/*
Initialisation of Formant_Frame is taken care of in Roots_into_Formant_Frame!
*/
try {
LPC_Frame_into_Formant_Frame (lpc, formant, my samplingPeriod, margin);
} catch (MelderError) {
Melder_clearError();
numberOfSuspectFrames ++;
}
if (interval == 1 || i % interval == 1)
Melder_progress ((double) i / my nx, U"LPC to Formant: frame ", i, U" out of ", my nx, U".");
if (interval == 1 || iframe % interval == 1)
Melder_progress ((double) iframe / my nx, U"LPC to Formant: frame ", iframe, U" out of ", my nx, U".");
}
Formant_sort (thee.get());
if (numberOfSuspectFrames > 0)
......@@ -121,10 +118,10 @@ autoFormant LPC_to_Formant (LPC me, double margin) {
}
void Formant_Frame_into_LPC_Frame (Formant_Frame me, LPC_Frame thee, double samplingPeriod) {
const double nyquistFrequency = 0.5 / samplingPeriod;
integer numberOfPoles = 2 * my nFormants;
if (my nFormants < 1)
return;
const double nyquistFrequency = 0.5 / samplingPeriod;
integer numberOfPoles = 2 * my nFormants;
autoVEC lpc = newVECzero (numberOfPoles + 2); // all odd coefficients have to be initialized to zero
lpc [2] = 1.0;
integer m = 2;
......@@ -160,8 +157,8 @@ autoLPC Formant_to_LPC (Formant me, double samplingPeriod) {
for (integer i = 1; i <= my nx; i ++) {
const Formant_Frame f = & my frames [i];
const LPC_Frame lpc = & thy d_frames [i];
const integer m = 2 * f -> nFormants; // TODO: what is m?
LPC_Frame_init (lpc, m);
const integer numberOfCoefficients = 2 * f -> nFormants;
LPC_Frame_init (lpc, numberOfCoefficients);
Formant_Frame_into_LPC_Frame (f, lpc, samplingPeriod);
}
return thee;
......
LPC/LPC_and_LFCC.cpp
View file @ 80504b54
/* LPC_and_LFCC.cpp
*
* Copyright (C) 1994-2018 David Weenink
* Copyright (C) 1994-2019 David Weenink
*
* This code is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -27,7 +27,8 @@
void LPC_Frame_into_CC_Frame (LPC_Frame me, CC_Frame thee) {
thy c0 = 0.5 * log (my gain);
if (my nCoefficients < 1) return;
if (my nCoefficients < 1)
return;
thy c [1] = - my a [1];
for (integer n = 2; n <= std::min ((integer) my nCoefficients, thy numberOfCoefficients); n ++) {
......@@ -45,17 +46,17 @@ void LPC_Frame_into_CC_Frame (LPC_Frame me, CC_Frame thee) {
}
void CC_Frame_into_LPC_Frame (CC_Frame me, LPC_Frame thee) {
integer n = std::min (my numberOfCoefficients, (integer) thy nCoefficients);
const integer numberOfCoefficients = std::min (my numberOfCoefficients, (integer) thy nCoefficients);
if (numberOfCoefficients < 1)
return;
thy gain = exp (2.0 * my c0);
if (n < 1) return;
thy a [1] = - my c [1];
for (integer i = 2; i <= n; i ++) {
for (integer i = 2; i <= numberOfCoefficients; i ++) {
longdouble ai = my c [i] * i;
for (integer j = 1; j < i; j ++)
ai += thy a [j] * my c [i - j] * (i - j);
thy a [i] = - double (ai) / i;
thy a [i] = - double (ai / i);
}
}
......@@ -66,9 +67,9 @@ autoLFCC LPC_to_LFCC (LPC me, integer numberOfCoefficients) {
autoLFCC thee = LFCC_create (my xmin, my xmax, my nx, my dx, my x1, numberOfCoefficients, 0, 0.5 / my samplingPeriod);
for (integer i = 1; i <= my nx; i ++) {
CC_Frame_init (& thy frame [i], numberOfCoefficients);
LPC_Frame_into_CC_Frame (& my d_frames [i], & thy frame [i]);
for (integer iframe = 1; iframe <= my nx; iframe ++) {
CC_Frame_init (& thy frame [iframe], numberOfCoefficients);
LPC_Frame_into_CC_Frame (& my d_frames [iframe], & thy frame [iframe]);
}
return thee;
} catch (MelderError) {
......@@ -84,9 +85,9 @@ autoLPC LFCC_to_LPC (LFCC me, integer numberOfCoefficients) {
numberOfCoefficients = std::min (numberOfCoefficients, my maximumNumberOfCoefficients);
autoLPC thee = LPC_create (my xmin, my xmax, my nx, my dx, my x1, numberOfCoefficients, 0.5 / my fmax);
for (integer i = 1; i <= my nx; i ++) {
LPC_Frame_init (& thy d_frames [i], numberOfCoefficients);
CC_Frame_into_LPC_Frame (& my frame [i], & thy d_frames [i]);
for (integer iframe = 1; iframe <= my nx; iframe ++) {
LPC_Frame_init (& thy d_frames [iframe], numberOfCoefficients);
CC_Frame_into_LPC_Frame (& my frame [iframe], & thy d_frames [iframe]);
}
return thee;
} catch (MelderError) {
......
LPC/LPC_and_LineSpectralFrequencies.cpp
View file @ 80504b54
......@@ -31,7 +31,6 @@
From: Joseph Rothweiler (1999), "On Polynomial Reduction in the Computation of LSP Frequencies."
IEEE Trans. on ASSP 7, 592--594.
*/
static void cos2x (VECVU const& g) {
for (integer i = 3; i <= g.size; i ++) {
for (integer j = g.size; j > i; j --)
......@@ -41,9 +40,10 @@ static void cos2x (VECVU const& g) {
}
static void Polynomial_fromLPC_Frame_lspsum (Polynomial me, LPC_Frame lpc) {
// Fs (z) = A(z) + z^-(p+1) A(1/z)
integer order = lpc -> nCoefficients, g_order = (order + 1) / 2; // orders
/*
Fs (z) = A(z) + z^-(p+1) A(1/z)
*/
const integer order = lpc -> nCoefficients, g_order = (order + 1) / 2; // orders
my coefficients [order + 2] = 1.0;
for (integer i = 1; i <= order; i ++)
my coefficients [order + 2 - i] = lpc -> a [i] + lpc -> a [order + 1 - i];
......@@ -53,20 +53,24 @@ static void Polynomial_fromLPC_Frame_lspsum (Polynomial me, LPC_Frame lpc) {
if (order % 2 == 0) // order even
Polynomial_divide_firstOrderFactor (me, -1.0, nullptr);
// transform to cos(w) terms
/*
Transform to cos(w) terms
*/
for (integer i = 1; i <= g_order + 1; i ++)
my coefficients [i] = my coefficients [g_order + i];
my numberOfCoefficients = g_order + 1;
// to Chebychev
/*
To Chebychev
*/
cos2x (my coefficients.part(1, my numberOfCoefficients));
}
static void Polynomial_fromLPC_Frame_lspdif (Polynomial me, LPC_Frame lpc) {
// Fa (z) = A(z) - z^-(p+1)A(1/z)
integer order = lpc -> nCoefficients;
/*
Fa (z) = A(z) - z^-(p+1)A(1/z)
*/
const integer order = lpc -> nCoefficients;
my coefficients [order + 2] = -1.0;
for (integer i = 1; i <= order; i ++)
my coefficients [order + 2 - i] = - lpc -> a [i] + lpc -> a [order + 1 - i];
......@@ -74,18 +78,28 @@ static void Polynomial_fromLPC_Frame_lspdif (Polynomial me, LPC_Frame lpc) {
my coefficients [1] = 1.0;
my numberOfCoefficients = order + 2;
if (order % 2 == 0) // Fa(z)/(z-1)
if (order % 2 == 0) {
/*
Fa(z)/(z-1)
*/
Polynomial_divide_firstOrderFactor (me, 1.0, nullptr);
else // Fa(z) / (z^2 - 1)
} else {
/*
Fa(z) / (z^2 - 1)
*/
Polynomial_divide_secondOrderFactor (me, 1.0);
// transform to cos(w) terms
}
/*
Transform to cos(w) terms
*/
integer g_order = my numberOfCoefficients / 2;
for (integer i = 1; i <= g_order + 1; i ++)
my coefficients [i] = my coefficients [g_order + i];
my numberOfCoefficients = g_order + 1;
// to Chebychev
/*
To Chebychev
*/
cos2x (my coefficients.part(1, my numberOfCoefficients));
}
......@@ -177,28 +191,27 @@ static integer Roots_fromPolynomial_grid (Roots me, Polynomial thee, double grid
while (xmin < thy xmax && numberOfRootsFound < thy numberOfCoefficients - 1) {
double xmax = xmin + gridSize;
xmax = xmax > thy xmax ? thy xmax : xmax;
//double root = Polynomial_findOneRealRoot_nr (thee, xmin, xmax);
double root = Polynomial_findOneSimpleRealRoot_ridders (thee, xmin, xmax);
const double root = Polynomial_findOneSimpleRealRoot_ridders (thee, xmin, xmax);
if (isdefined (root) && (numberOfRootsFound == 0 || my v [numberOfRootsFound].re != root)) {
my v [++ numberOfRootsFound].re = root; // root not at border of interval
my v [numberOfRootsFound].im = 0.0;
}
xmin = xmax;
}
// make rest of storage undefined (not necessary)
return numberOfRootsFound;
}
static void LineSpectralFrequencies_Frame_initFromLPC_Frame_grid (LineSpectralFrequencies_Frame me, LPC_Frame thee, Polynomial g1, Polynomial g2, Roots roots, double gridSize, double maximumFrequency) {
/* Construct Fs and Fa
/*
Construct Fs and Fa
divide out the zeros
transform to polynomial equations g1 and g2 of half the order
*/
LineSpectralFrequencies_Frame_init (me, thy nCoefficients);
Polynomial_fromLPC_Frame_lspsum (g1, thee);
integer half_order_g1 = g1 -> numberOfCoefficients - 1;
const integer half_order_g1 = g1 -> numberOfCoefficients - 1;
Polynomial_fromLPC_Frame_lspdif (g2, thee);
integer half_order_g2 = g2 -> numberOfCoefficients - 1;
const integer half_order_g2 = g2 -> numberOfCoefficients - 1;
integer numberOfBisections = 0, numberOfRootsFound = 0;
while (numberOfRootsFound < half_order_g1 && numberOfBisections < 10) {
......@@ -207,19 +220,21 @@ static void LineSpectralFrequencies_Frame_initFromLPC_Frame_grid (LineSpectralFr
numberOfBisections++;
}
Melder_require (numberOfBisections < 10, U"Too many bisections.");
// [g1-> xmin, g1 -> xmax] <==> [nyquistFrequency, 0] i.e. highest root corresponds to lowest frequency
Melder_require (numberOfBisections < 10,
U"Too many bisections.");
/*
[g1-> xmin, g1 -> xmax] <==> [nyquistFrequency, 0],
i.e. highest root corresponds to lowest frequency
*/
for (integer i = 1; i <= half_order_g1; i ++)
my frequencies [2 * i - 1] = acos (roots -> v [half_order_g1 + 1 - i].re / 2.0) / NUMpi * maximumFrequency;
// the roots of g2 lie inbetween the roots of g1
/*
The roots of g2 lie inbetween the roots of g1
*/
for (integer i = 1; i <= half_order_g2; i ++) {
double xmax = roots -> v [half_order_g1 + 1 - i].re;
double xmin = i == half_order_g1 ? g1 -> xmin : roots -> v [half_order_g1 - i].re;
double root = Polynomial_findOneSimpleRealRoot_ridders (g2, xmin, xmax);
const double xmax = roots -> v [half_order_g1 + 1 - i].re;
const double xmin = i == half_order_g1 ? g1 -> xmin : roots -> v [half_order_g1 - i].re;
const double root = Polynomial_findOneSimpleRealRoot_ridders (g2, xmin, xmax);
if (isdefined (root))
my frequencies [2 * i] = acos (root / 2.0) / NUMpi * maximumFrequency;
else
......@@ -232,16 +247,17 @@ autoLineSpectralFrequencies LPC_to_LineSpectralFrequencies (LPC me, double gridS
if (gridSize == 0.0)
gridSize = 0.02;
double nyquistFrequency = 0.5 / my samplingPeriod;
const double nyquistFrequency = 0.5 / my samplingPeriod;
autoLineSpectralFrequencies thee = LineSpectralFrequencies_create (my xmin, my xmax, my nx, my dx, my x1, my maxnCoefficients, nyquistFrequency);
autoPolynomial g1 = Polynomial_create (-2.0, 2.0, my maxnCoefficients + 1); // large enough
autoPolynomial g2 = Polynomial_create (-2.0, 2.0, my maxnCoefficients + 1);
autoRoots roots = Roots_create ((my maxnCoefficients + 1) / 2);
for (integer iframe = 1; iframe <= my nx; iframe ++) {
LPC_Frame lpc_frame = & my d_frames [iframe];
LineSpectralFrequencies_Frame lsf_frame = & thy d_frames [iframe];
/* Construct Fs and Fa
const LPC_Frame lpc_frame = & my d_frames [iframe];
const LineSpectralFrequencies_Frame lsf_frame = & thy d_frames [iframe];
/*
Construct Fs and Fa
divide out the zeros
transform to polynomial equations g1 and g2 of half the order
find zeros
......@@ -259,24 +275,22 @@ autoLineSpectralFrequencies LPC_to_LineSpectralFrequencies (LPC me, double gridS
*/
static void LPC_Frame_initFromLineSpectralFrequencies_Frame (LPC_Frame me, LineSpectralFrequencies_Frame thee, Polynomial fs, Polynomial fa, double maximumFrequency) {
LPC_Frame_init (me, thy numberOfFrequencies);
/*
Reconstruct Fs (z)
Use my a as a buffer whose size changes!!!
*/
integer numberOfOmegas = (thy numberOfFrequencies + 1) / 2;
for (integer i = 1; i <= numberOfOmegas; i ++) {
double omega = thy frequencies [2 * i -1] / maximumFrequency * NUMpi;
const double omega = thy frequencies [2 * i -1] / maximumFrequency * NUMpi;
my a [i] = -2.0 * cos (omega);
}
Polynomial_initFromProductOfSecondOrderTerms (fs, my a.part (1, numberOfOmegas));
/*
Reconstruct Fa (z)
*/
numberOfOmegas = thy numberOfFrequencies / 2;
for (integer i = 1; i <= numberOfOmegas; i ++) {
double omega = thy frequencies [2 * i] / maximumFrequency * NUMpi;
const double omega = thy frequencies [2 * i] / maximumFrequency * NUMpi;
my a [i] = -2.0 * cos (omega);
}
Polynomial_initFromProductOfSecondOrderTerms (fa, my a.part (1, numberOfOmegas));
......@@ -302,9 +316,10 @@ autoLPC LineSpectralFrequencies_to_LPC (LineSpectralFrequencies me) {
autoPolynomial fa = Polynomial_create (-1.0, 1.0, my maximumNumberOfFrequencies + 2);
for (integer iframe = 1; iframe <= my nx; iframe ++) {
LineSpectralFrequencies_Frame lsf_frame = & my d_frames [iframe];
LPC_Frame lpc_frame = & thy d_frames [iframe];
/* Construct Fs and Fa
const LineSpectralFrequencies_Frame lsf_frame = & my d_frames [iframe];
const LPC_Frame lpc_frame = & thy d_frames [iframe];
/*
Construct Fs and Fa
A(z) = (Fs(z) + Fa(z))/2
*/
LPC_Frame_initFromLineSpectralFrequencies_Frame (lpc_frame, lsf_frame, fs.get(), fa.get(), my maximumFrequency);
......
LPC/LPC_and_Polynomial.cpp
View file @ 80504b54
/* LPC_and_Polynomial.cpp
*
* Copyright (C) 1994-2018 David Weenink
* Copyright (C) 1994-2019 David Weenink
*
* This code is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -23,7 +23,7 @@
#include "LPC_and_Polynomial.h"
autoPolynomial LPC_Frame_to_Polynomial (LPC_Frame me) {
integer degree = my nCoefficients;
const integer degree = my nCoefficients;
autoPolynomial thee = Polynomial_create (-1, 1, degree);
for (integer i = 1; i <= degree; i ++)
thy coefficients [i] = my a [degree - i + 1];
......@@ -34,7 +34,7 @@ autoPolynomial LPC_Frame_to_Polynomial (LPC_Frame me) {
autoPolynomial LPC_to_Polynomial (LPC me, double time) {
try {
integer iFrame = Sampled_xToIndex (me, time);
Melder_clip (integer (1), & iFrame, my nx); // constant extrapolation
Melder_clip (1_integer, & iFrame, my nx); // constant extrapolation
autoPolynomial thee = LPC_Frame_to_Polynomial (& my d_frames [iFrame]);
return thee;
} catch (MelderError) {
......
LPC/LPC_and_Tube.cpp
View file @ 80504b54
......@@ -43,7 +43,9 @@ void LPC_Frame_into_Tube_Frame_area (LPC_Frame me, Tube_Frame thee) {
double VocalTract_LPC_Frame_getMatchingLength (VocalTract me, LPC_Frame thee, double glottalDamping, bool radiationDamping, bool internalDamping) {
try {
// match the average distance between the first two formants in the VocaTract and the LPC spectrum
/*
Match the average distance between the first two formants in the VocaTract and the LPC spectrum
*/
const integer numberOfFrequencies = 1000;
const double maximumFrequency = 5000.0;
autoSpectrum vts = VocalTract_to_Spectrum (me, numberOfFrequencies, maximumFrequency, glottalDamping, radiationDamping, internalDamping);
......@@ -70,8 +72,9 @@ double LPC_Frame_getVTL_wakita (LPC_Frame me, double samplingPeriod, double refL
struct structTube_Frame rc_struct, af_struct;
const Tube_Frame rc = & rc_struct, af = & af_struct;
try {
const integer m = my nCoefficients;
double length, dlength = 0.001, wakita_length = undefined;
const integer numberOfCoefficients = my nCoefficients;
const double dlength = 0.001;
double length, wakita_length = undefined;
double varMin = 1e308;
memset (& lpc_struct, 0, sizeof (lpc_struct));
......@@ -80,31 +83,30 @@ double LPC_Frame_getVTL_wakita (LPC_Frame me, double samplingPeriod, double refL
memset (& af_struct, 0, sizeof (af_struct));
LPC_Frame_init (lpc, m);
Tube_Frame_init (rc, m, refLength);
Tube_Frame_init (af, m, refLength);
// Step 2
LPC_Frame_init (lpc, numberOfCoefficients);
Tube_Frame_init (rc, numberOfCoefficients, refLength);
Tube_Frame_init (af, numberOfCoefficients, refLength);
/*
Step 2
*/
LPC_Frame_into_Formant_Frame (me, f, samplingPeriod, 0);
// LPC_Frame_into_Formant_Frame performs the Formant_Frame_init !!
Melder_require (f -> nFormants > 0, U"Not enough formants.");
/*
LPC_Frame_into_Formant_Frame performs the Formant_Frame_init !!
*/
Melder_require (f -> nFormants > 0,
U"Not enough formants.");
double *area = af -> c.at; // TODO
double lmin = length = 0.10;
double plength = refLength;
while (length <= 0.25) {
// Step 3
/*
Step 3
*/
const double fscale = plength / length;
for (integer i = 1; i <= f -> nFormants; i ++) {
f -> formant [i]. frequency *= fscale;
f -> formant [i]. bandwidth *= fscale;
}
/*
20000125: Bovenstaande schaling van f1/b1 kan ook gedaan worden door
MGfb_to_a (f, b, nf, samplingFrequency*length/refLength, a1)
......@@ -115,35 +117,35 @@ double LPC_Frame_getVTL_wakita (LPC_Frame me, double samplingPeriod, double refL
refLength=c*aantalFormanten/Fs waarbij c=340 m/s (geluidssnelheid).
Bij Fs=10000 zou aantalFormanten=5 zijn en refLength -> 0.17 m
*/
// step 4
/*
Step 4
*/
Formant_Frame_into_LPC_Frame (f, lpc, samplingPeriod);
// step 5
/*
Step 5
*/
rc -> length = length;
LPC_Frame_into_Tube_Frame_rc (lpc, rc);
// step 6.1
/*
Step 6.1
*/
Tube_Frames_rc_into_area (rc, af);
// step 6.2 Log(areas)
/*
Step 6.2 Log(areas)
*/
double logSum = 0.0;
for (integer i = 1; i <= af -> numberOfSegments; i ++) {
area [i] = log (area [i]);
logSum += area [i];
}
// step 6.3 and 7
/*
Steps 6.3 and 7
*/
double var = 0.0;
for (integer i = 1; i <= af -> numberOfSegments; i ++) {
const double delta = area [i] - logSum / af -> numberOfSegments;
var += delta * delta;
}
if (var < varMin) {
lmin = length;
varMin = var;
......@@ -187,8 +189,8 @@ void VocalTract_setLength (VocalTract me, double newLength) {
autoVocalTract LPC_to_VocalTract_slice_special (LPC me, double time, double glottalDamping, bool radiationDamping, bool internalDamping) {
try {
integer frameNumber = Sampled_xToLowIndex (me, time); // ppgb: BUG? Is rounding down the correct thing to do? not nearestIndex?
Melder_clip (integer (1), & frameNumber, my nx); // constant extrapolation
integer frameNumber = Sampled_xToNearestIndex (me, time);
Melder_clip (1_integer, & frameNumber, my nx); // constant extrapolation
LPC_Frame lpc = & my d_frames [frameNumber];
autoVocalTract thee = LPC_Frame_to_VocalTract (lpc, 0.17);
const double length = VocalTract_LPC_Frame_getMatchingLength (thee.get(), lpc, glottalDamping, radiationDamping, internalDamping);
......@@ -217,8 +219,8 @@ autoVocalTract LPC_Frame_to_VocalTract (LPC_Frame me, double length) {
autoVocalTract LPC_to_VocalTract_slice (LPC me, double time, double length) {
try {
integer frameNumber = Sampled_xToNearestIndex (me, time);
Melder_clip (integer (1), & frameNumber, my nx); // constant extrapolation
LPC_Frame lpc = & my d_frames [frameNumber];
Melder_clip (1_integer, & frameNumber, my nx); // constant extrapolation
const LPC_Frame lpc = & my d_frames [frameNumber];
autoVocalTract thee = LPC_Frame_to_VocalTract (lpc, length);
return thee;
} catch (MelderError) {
......
LPC/LPC_to_Spectrogram.cpp
View file @ 80504b54
......@@ -37,12 +37,12 @@ autoSpectrogram LPC_to_Spectrogram (LPC me, double dfMin, double bandwidthReduct
autoSpectrogram thee = Spectrogram_create (my xmin, my xmax, my nx, my dx, my x1, 0.0, samplingFrequency / 2.0, nfft / 2 + 1, freqStep, 0.0);
for (integer i = 1; i <= my nx; i ++) {
const double t = Sampled_indexToX (me, i);
for (integer iframe = 1; iframe <= my nx; iframe ++) {
const double t = Sampled_indexToX (me, iframe);
autoSpectrum spec = LPC_to_Spectrum (me, t, dfMin, bandwidthReduction, deEmphasisFrequency);
for (integer j = 1; j <= spec -> nx; j ++) {
const double re = spec -> z [1] [j], im = spec -> z [2] [j];
thy z [j] [i] = re * re + im * im;
thy z [j] [iframe] = re * re + im * im;
}
}
return thee;
......