Skip to content
Commits on Source (11)
Expand all Show whitespace changes
Inline Side-by-side
EEG/EEG.cpp
View file @ afb7171d
......@@ -409,18 +409,16 @@ autoEEG EEG_readFromBdfFile (MelderFile file) {
}
}
static void detrend (double *a, integer numberOfSamples) {
double firstValue = a [1], lastValue = a [numberOfSamples];
a [1] = a [numberOfSamples] = 0.0;
for (integer isamp = 2; isamp < numberOfSamples; isamp ++) {
a [isamp] -= ((isamp - 1.0) * lastValue + (numberOfSamples - isamp) * firstValue) / (numberOfSamples - 1);
}
static void detrend (VEC const& channel) {
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);
}
void EEG_detrend (EEG me) {
for (integer ichan = 1; ichan <= my numberOfChannels - EEG_getNumberOfExtraSensors (me); ichan ++) {
detrend (my sound -> z [ichan], my sound -> nx);
}
for (integer ichan = 1; ichan <= my numberOfChannels - EEG_getNumberOfExtraSensors (me); ichan ++)
detrend (my sound -> z.row (ichan));
}
void EEG_filter (EEG me, double lowFrequency, double lowWidth, double highFrequency, double highWidth, bool doNotch50Hz) {
......@@ -441,7 +439,7 @@ void EEG_filter (EEG me, double lowFrequency, double lowWidth, double highFreque
Spectrum_stopHannBand (spec.get(), 48.0, 52.0, 1.0);
}
autoSound him = Spectrum_to_Sound (spec.get());
NUMvector_copyElements (his z [1], my sound -> z [ichan], 1, my sound -> nx);
NUMvector_copyElements (& his z [1] [0], & my sound -> z [ichan] [0], 1, my sound -> nx);
}
} catch (MelderError) {
Melder_throw (me, U": not filtered.");
......@@ -496,25 +494,22 @@ void EEG_subtractMeanChannel (EEG me, integer fromChannel, integer toChannel) {
const integer numberOfElectrodeChannels = my numberOfChannels - EEG_getNumberOfExtraSensors (me);
for (integer isamp = 1; isamp <= my sound -> nx; isamp ++) {
double referenceValue = 0.0;
for (integer ichan = fromChannel; ichan <= toChannel; ichan ++) {
for (integer ichan = fromChannel; ichan <= toChannel; ichan ++)
referenceValue += my sound -> z [ichan] [isamp];
}
referenceValue /= (toChannel - fromChannel + 1);
for (integer ichan = 1; ichan <= numberOfElectrodeChannels; ichan ++) {
for (integer ichan = 1; ichan <= numberOfElectrodeChannels; ichan ++)
my sound -> z [ichan] [isamp] -= referenceValue;
}
}
}
void EEG_setChannelToZero (EEG me, integer channelNumber) {
try {
if (channelNumber < 1 || channelNumber > my numberOfChannels)
Melder_throw (U"No channel ", channelNumber, U".");
integer numberOfSamples = my sound -> nx;
double *channel = my sound -> z [channelNumber];
for (integer isample = 1; isample <= numberOfSamples; isample ++) {
VEC channel = my sound -> z.row (channelNumber);
for (integer isample = 1; isample <= numberOfSamples; isample ++)
channel [isample] = 0.0;
}
} catch (MelderError) {
Melder_throw (me, U": channel ", channelNumber, U" not set to zero.");
}
......@@ -594,9 +589,8 @@ static void Sound_removeChannel (Sound me, integer channelNumber) {
U"No channel ", channelNumber, U".");
Melder_require (my ny > 1,
U"Cannot remove last remaining channel.");
for (integer ichan = channelNumber; ichan < my ny; ichan ++) {
NUMvector_copyElements (my z [ichan + 1], my z [ichan], 1, my nx);
}
for (integer ichan = channelNumber; ichan < my ny; ichan ++)
NUMvector_copyElements (& my z [ichan + 1] [0], & my z [ichan] [0], 1, my nx);
my ymax -= 1.0;
my ny -= 1;
} catch (MelderError) {
......@@ -636,7 +630,7 @@ autoEEG EEGs_concatenate (OrderedOf<structEEG>* me) {
Melder_throw (U"Cannot concatenate zero EEG objects.");
EEG first = my at [1];
integer numberOfChannels = first -> numberOfChannels;
autostring32vector channelNames = STRVECclone (first -> channelNames.get());
autostring32vector channelNames = newSTRVECcopy (first -> channelNames.get());
for (integer ieeg = 2; ieeg <= my size; ieeg ++) {
EEG other = my at [ieeg];
if (other -> numberOfChannels != numberOfChannels)
......@@ -671,7 +665,7 @@ autoEEG EEG_extractPart (EEG me, double tmin, double tmax, bool preserveTimes) {
try {
autoEEG thee = Thing_new (EEG);
thy numberOfChannels = my numberOfChannels;
thy channelNames = STRVECclone (my channelNames.get());
thy channelNames = newSTRVECcopy (my channelNames.get());
thy sound = Sound_extractPart (my sound.get(), tmin, tmax, kSound_windowShape::RECTANGULAR, 1.0, preserveTimes);
thy textgrid = TextGrid_extractPart (my textgrid.get(), tmin, tmax, preserveTimes);
thy xmin = thy textgrid -> xmin;
......@@ -723,7 +717,7 @@ autoEEG EEG_MixingMatrix_to_EEG_unmix (EEG me, MixingMatrix you) {
his sound = Sound_MixingMatrix_unmix (my sound.get(), you);
his textgrid = Data_copy (my textgrid.get());
his numberOfChannels = your numberOfColumns;
his channelNames = STRVECclone (your columnLabels.get());
his channelNames = newSTRVECcopy (your columnLabels.get());
return him;
}
......@@ -741,7 +735,7 @@ autoEEG EEG_MixingMatrix_to_EEG_mix (EEG me, MixingMatrix you) {
his sound = Sound_MixingMatrix_mix (my sound.get(), you);
his textgrid = Data_copy (my textgrid.get());
his numberOfChannels = your numberOfRows;
his channelNames = STRVECclone (your rowLabels.get());
his channelNames = newSTRVECcopy (your rowLabels.get());
return him;
}
......
EEG/ERP.cpp
View file @ afb7171d
......@@ -78,7 +78,7 @@ void ERP_drawChannel_number (ERP me, Graphics graphics, integer channelNumber, d
*/
Graphics_setInner (graphics);
Graphics_setWindow (graphics, tmin, tmax, vmin, vmax);
Graphics_function (graphics, my z [channelNumber], ixmin, ixmax, Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax));
Graphics_function (graphics, & my z [channelNumber] [0], ixmin, ixmax, Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax));
Graphics_unsetInner (graphics);
if (garnish) {
Graphics_drawInnerBox (graphics);
......
EEG/ERPTier.cpp
View file @ afb7171d
......@@ -70,7 +70,7 @@ static autoERPTier EEG_PointProcess_to_ERPTier (EEG me, PointProcess events, dou
Function_init (thee.get(), fromTime, toTime);
thy numberOfChannels = my numberOfChannels - EEG_getNumberOfExtraSensors (me);
Melder_assert (thy numberOfChannels > 0);
thy channelNames = STRVECclone (my channelNames.get());
thy channelNames = newSTRVECcopy (my channelNames.get());
integer numberOfEvents = events -> nt;
double soundDuration = toTime - fromTime;
double samplingPeriod = my sound -> dx;
......@@ -191,7 +191,7 @@ void ERPTier_subtractBaseline (ERPTier me, double tmin, double tmax) {
ERPPoint event = my points.at [ievent];
for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
double mean = Vector_getMean (event -> erp.get(), tmin, tmax, ichannel);
double *channel = event -> erp -> z [ichannel];
VEC channel = event -> erp -> z.row (ichannel);
for (integer isample = 1; isample <= numberOfSamples; isample ++) {
channel [isample] -= mean;
}
......@@ -213,7 +213,7 @@ void ERPTier_rejectArtefacts (ERPTier me, double threshold) {
double minimum = event -> erp -> z [1] [1];
double maximum = minimum;
for (integer ichannel = 1; ichannel <= (numberOfChannels & ~ 15); ichannel ++) {
double *channel = event -> erp -> z [ichannel];
constVEC channel = event -> erp -> z.row (ichannel);
for (integer isample = 1; isample <= numberOfSamples; isample ++) {
double value = channel [isample];
if (value < minimum) minimum = value;
......@@ -233,19 +233,10 @@ autoERP ERPTier_extractERP (ERPTier me, integer eventNumber) {
Melder_throw (U"No events.");
ERPTier_checkEventNumber (me, eventNumber);
ERPPoint event = my points.at [eventNumber];
integer numberOfChannels = event -> erp -> ny;
Melder_assert (numberOfChannels == my numberOfChannels);
integer numberOfSamples = event -> erp -> nx;
Melder_assert (event -> erp -> ny == my numberOfChannels);
autoERP thee = Thing_new (ERP);
event -> erp -> structSound :: v_copy (thee.get());
for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
double *oldChannel = event -> erp -> z [ichannel];
double *newChannel = thy z [ichannel];
for (integer isample = 1; isample <= numberOfSamples; isample ++) {
newChannel [isample] = oldChannel [isample];
}
}
thy channelNames = STRVECclone (my channelNames.get());
thy channelNames = newSTRVECcopy (my channelNames.get());
return thee;
} catch (MelderError) {
Melder_throw (me, U": ERP not extracted.");
......@@ -258,29 +249,17 @@ autoERP ERPTier_to_ERP_mean (ERPTier me) {
if (numberOfEvents < 1)
Melder_throw (U"No events.");
ERPPoint firstEvent = my points.at [1];
integer numberOfChannels = firstEvent -> erp -> ny;
integer numberOfSamples = firstEvent -> erp -> nx;
Melder_assert (firstEvent -> erp -> ny == my numberOfChannels);
autoERP mean = Thing_new (ERP);
firstEvent -> erp -> structSound :: v_copy (mean.get());
Melder_assert (mean -> ny == my numberOfChannels);
for (integer ievent = 2; ievent <= numberOfEvents; ievent ++) {
ERPPoint event = my points.at [ievent];
for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
double *erpChannel = event -> erp -> z [ichannel];
double *meanChannel = mean -> z [ichannel];
for (integer isample = 1; isample <= numberOfSamples; isample ++) {
meanChannel [isample] += erpChannel [isample];
}
}
}
double factor = 1.0 / numberOfEvents;
for (integer ichannel = 1; ichannel <= numberOfChannels; ichannel ++) {
double *meanChannel = mean -> z [ichannel];
for (integer isample = 1; isample <= numberOfSamples; isample ++) {
meanChannel [isample] *= factor;
Melder_assert (event -> erp -> ny == my numberOfChannels);
mean -> z.all() += event -> erp -> z.all();
}
}
Melder_assert (mean -> ny == my numberOfChannels);
mean -> channelNames = STRVECclone (my channelNames.get());
mean -> z.all() *= 1.0 / numberOfEvents;
mean -> channelNames = newSTRVECcopy (my channelNames.get());
return mean;
} catch (MelderError) {
Melder_throw (me, U": mean not computed.");
......@@ -297,7 +276,7 @@ autoERPTier ERPTier_extractEventsWhereColumn_number (ERPTier me, Table table, in
autoERPTier thee = Thing_new (ERPTier);
Function_init (thee.get(), my xmin, my xmax);
thy numberOfChannels = my numberOfChannels;
thy channelNames = STRVECclone (my channelNames.get());
thy channelNames = newSTRVECcopy (my channelNames.get());
for (integer ievent = 1; ievent <= my points.size; ievent ++) {
ERPPoint oldEvent = my points.at [ievent];
TableRow row = table -> rows.at [ievent];
......@@ -326,7 +305,7 @@ autoERPTier ERPTier_extractEventsWhereColumn_string (ERPTier me, Table table,
autoERPTier thee = Thing_new (ERPTier);
Function_init (thee.get(), my xmin, my xmax);
thy numberOfChannels = my numberOfChannels;
thy channelNames = STRVECclone (my channelNames.get());
thy channelNames = newSTRVECcopy (my channelNames.get());
for (integer ievent = 1; ievent <= my points.size; ievent ++) {
ERPPoint oldEvent = my points.at [ievent];
TableRow row = table -> rows.at [ievent];
......
EEG/ERPWindow.cpp
View file @ afb7171d
......@@ -158,7 +158,7 @@ void ERP_drawScalp_garnish (Graphics graphics, double vmin, double vmax, enum kG
}
}
Graphics_setColourScale (graphics, colourScale);
Graphics_image (graphics, legend.at, 1, 2, 0.85, 0.98, 1, n, -0.8, +0.8, 0.0, 1.0);
Graphics_image (graphics, legend.all(), 0.85, 0.98, -0.8, +0.8, 0.0, 1.0);
Graphics_setColourScale (graphics, kGraphics_colourScale::GREY);
Graphics_rectangle (graphics, 0.85, 0.98, -0.8, +0.8);
Graphics_setTextAlignment (graphics, Graphics_RIGHT, Graphics_TOP);
......@@ -232,7 +232,7 @@ void ERP_drawScalp (ERP me, Graphics graphics, double tmin, double tmax, double
}
}
}
Graphics_image (graphics, image.at, 1, n, -1.0-0.5/n, 1.0+0.5/n, 1, n, -1.0-0.5/n, 1.0+0.5/n, vmin, vmax);
Graphics_image (graphics, image.all(), -1.0-0.5/n, 1.0+0.5/n, -1.0-0.5/n, 1.0+0.5/n, vmin, vmax);
Graphics_setColourScale (graphics, kGraphics_colourScale::GREY);
Graphics_setLineWidth (graphics, 2.0);
/*
......@@ -357,7 +357,7 @@ void structERPWindow :: v_drawSelectionViewer () {
}
}
Graphics_setColourScale (our graphics.get(), our p_scalp_colourScale);
Graphics_image (our graphics.get(), image.at, 1, n, -1.0-0.5/n, 1.0+0.5/n, 1, n, -1.0-0.5/n, 1.0+0.5/n, minimum, maximum);
Graphics_image (our graphics.get(), image.all(), -1.0-0.5/n, 1.0+0.5/n, -1.0-0.5/n, 1.0+0.5/n, minimum, maximum);
Graphics_setColourScale (our graphics.get(), kGraphics_colourScale::GREY);
Graphics_setLineWidth (our graphics.get(), 2.0);
/*
......@@ -393,7 +393,8 @@ void structERPWindow :: v_drawSelectionViewer () {
OPTIONMENU_ENUM_VARIABLE (kGraphics_colourScale, v_prefs_scalpColourSpace)
void structERPWindow :: v_prefs_addFields (EditorCommand cmd) {
OPTIONMENU_ENUM_FIELD (v_prefs_scalpColourSpace, U"Scalp colour space", kGraphics_colourScale, kGraphics_colourScale::BLUE_TO_RED)
OPTIONMENU_ENUM_FIELD (kGraphics_colourScale, v_prefs_scalpColourSpace,
U"Scalp colour space", kGraphics_colourScale::BLUE_TO_RED)
}
void structERPWindow :: v_prefs_setValues (EditorCommand cmd) {
SET_ENUM (v_prefs_scalpColourSpace, kGraphics_colourScale, p_scalp_colourScale)
......
EEG/praat_EEG.cpp
View file @ afb7171d
......@@ -125,12 +125,13 @@ DO
}
FORM (MODIFY_EEG_removeTriggers, U"Remove triggers", nullptr) {
OPTIONMENU_ENUM (removeEveryTriggerThat___, U"Remove every trigger that...", kMelder_string, DEFAULT)
OPTIONMENU_ENUM (kMelder_string, removeEveryTriggerThat___,
U"Remove every trigger that...", kMelder_string::DEFAULT)
SENTENCE (___theText, U"...the text", U"hi")
OK
DO
MODIFY_EACH (EEG)
EEG_removeTriggers (me, (kMelder_string) removeEveryTriggerThat___, ___theText);
EEG_removeTriggers (me, removeEveryTriggerThat___, ___theText);
MODIFY_EACH_END
}
......@@ -266,27 +267,30 @@ DO
FORM (NEW_EEG_to_ERPTier_triggers, U"To ERPTier (triggers)", nullptr) {
REAL (fromTime, U"From time (s)", U"-0.11")
REAL (toTime, U"To time (s)", U"0.39")
OPTIONMENU_ENUM (getEveryEventWithATriggerThat, U"Get every event with a trigger that", kMelder_string, DEFAULT)
OPTIONMENU_ENUM (kMelder_string, getEveryEventWithATriggerThat,
U"Get every event with a trigger that", kMelder_string::DEFAULT)
SENTENCE (theText, U"...the text", U"1")
OK
DO
CONVERT_EACH (EEG)
autoERPTier result = EEG_to_ERPTier_triggers (me, fromTime, toTime, (kMelder_string) getEveryEventWithATriggerThat, theText);
autoERPTier result = EEG_to_ERPTier_triggers (me, fromTime, toTime, getEveryEventWithATriggerThat, theText);
CONVERT_EACH_END (my name.get(), U"_trigger", theText)
}
FORM (NEW_EEG_to_ERPTier_triggers_preceded, U"To ERPTier (triggers, preceded)", nullptr) {
REAL (fromTime, U"From time (s)", U"-0.11")
REAL (toTime, U"To time (s)", U"0.39")
OPTIONMENU_ENUM (getEveryEventWithATriggerThat, U"Get every event with a trigger that", kMelder_string, DEFAULT)
OPTIONMENU_ENUM (kMelder_string, getEveryEventWithATriggerThat,
U"Get every event with a trigger that", kMelder_string::DEFAULT)
SENTENCE (text1, U"...the text", U"1")
OPTIONMENU_ENUM (andIsPrecededByATriggerThat, U"and is preceded by a trigger that", kMelder_string, DEFAULT)
OPTIONMENU_ENUM (kMelder_string, andIsPrecededByATriggerThat,
U"and is preceded by a trigger that", kMelder_string::DEFAULT)
SENTENCE (text2, U" ...the text", U"4")
OK
DO
CONVERT_EACH (EEG)
autoERPTier result = EEG_to_ERPTier_triggers_preceded (me, fromTime, toTime,
(kMelder_string) getEveryEventWithATriggerThat, text1, (kMelder_string) andIsPrecededByATriggerThat, text2);
(kMelder_string) getEveryEventWithATriggerThat, text1, andIsPrecededByATriggerThat, text2);
CONVERT_EACH_END (my name.get(), U"_trigger", text2)
}
......@@ -422,7 +426,8 @@ FORM (GRAPHICS_ERP_drawScalp_colour, U"ERP: Draw scalp (colour)", nullptr) {
REAL (toTime, U"right Time range", U"0.2")
REAL (fromVoltage, U"left Voltage range (V)", U"10e-6")
REAL (toVoltage, U"right Voltage range", U"-10e-6")
RADIO_ENUM (colourScale, U"Colour scale", kGraphics_colourScale, BLUE_TO_RED)
RADIO_ENUM (kGraphics_colourScale, colourScale,
U"Colour scale", kGraphics_colourScale::BLUE_TO_RED)
BOOLEAN (garnish, U"Garnish", true)
OK
DO
......@@ -435,7 +440,8 @@ DO
FORM (GRAPHICS_ERP_drawScalp_garnish, U"ERP: Draw scalp (garnish)", nullptr) {
REAL (fromVoltage, U"left Voltage range (V)", U"10e-6")
REAL (toVoltage, U"right Voltage range", U"-10e-6")
RADIO_ENUM (colourScale, U"Colour scale", kGraphics_colourScale, BLUE_TO_RED)
RADIO_ENUM (kGraphics_colourScale, colourScale,
U"Colour scale", kGraphics_colourScale::BLUE_TO_RED)
OK
DO
GRAPHICS_NONE
......@@ -704,7 +710,7 @@ DIRECT (NEW_ERPTier_to_ERP_mean) {
FORM (NEW1_ERPTier_Table_extractEventsWhereColumn_number, U"Extract events where column (number)", nullptr) {
WORD (extractAllEventsWhereColumn___, U"Extract all events where column...", U"")
RADIO_ENUM (___is___, U"...is...", kMelder_number, DEFAULT)
RADIO_ENUM (kMelder_number, ___is___, U"...is...", kMelder_number::DEFAULT)
REAL (___theNumber, U"...the number", U"0.0")
OK
DO
......@@ -716,13 +722,13 @@ DO
FORM (NEW1_ERPTier_Table_extractEventsWhereColumn_text, U"Extract events where column (text)", nullptr) {
WORD (extractAllEventsWhereColumn___, U"Extract all events where column...", U"")
OPTIONMENU_ENUM (___, U"...", kMelder_string, DEFAULT)
OPTIONMENU_ENUM (kMelder_string, ___, U"...", kMelder_string::DEFAULT)
SENTENCE (___theText, U"...the text", U"hi")
OK
DO
CONVERT_TWO (ERPTier, Table)
integer columnNumber = Table_getColumnIndexFromColumnLabel (you, extractAllEventsWhereColumn___);
autoERPTier result = ERPTier_extractEventsWhereColumn_string (me, you, columnNumber, (kMelder_string) ___, ___theText);
autoERPTier result = ERPTier_extractEventsWhereColumn_string (me, you, columnNumber, ___, ___theText);
CONVERT_TWO_END (my name.get())
}
......
FFNet/FFNet.cpp
View file @ afb7171d
This diff is collapsed.
FFNet/FFNet.h
View file @ afb7171d
......@@ -2,7 +2,7 @@
#define _FFNet_h_
/* FFNet.h
*
* Copyright (C) 1997-207 David Weenink
* Copyright (C) 1997-2018 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
......@@ -34,48 +34,48 @@
#include "Minimizers.h"
#include "PatternList.h"
#include "TableOfReal.h"
#include "melder.h"
#include "FFNet_def.h"
/* Parameters:
* integer nLayers : the #layers in the net (exclusive the inputs)
* nUnitsInLayer : array[0..nLayers] the #units in each layer
* nUnitsInLayer[0] : #inputs
* nUnitsInLayer[nLayers] :#outputs
* nWeights : the total #weights in the net (inclusive bias)
* double *w : array[1..nWeights] with the connection strengths
* *activity : array[1..nNodes] with activities
* integer numberOfLayers : the #layers in the net (exclusive the inputs)
* numberOfUnitsInLayer : array[1..numberOfLayers] the #units in each layer
* numberOfUnitsInLayer[numberOfLayers] :#outputs
* numberOfWeights : the total #weights in the net (inclusive bias)
* double *w : array[1..numberOfWeights] with the connection strengths
* *activity : array[1..numberOfNodes] with activities
* outputLabels : labels belonging to the outputs
* BOOKKEEPING:
* integer nNodes : total #nodes: bias modelled as unit with constant activity)
* *isbias : array[1..nNodes] set 1 if node is bias else 0
* *nodeFirst : array[1..nNodes] first node connected to this unit
* *nodeLast: : array[1..nNodes] last node connected to this unit
* *wFirst : array[1..nNodes] first index in *w for this unit
* *wLast : array[1..nNodes] last (inclusive the bias)
* integer numberOfNodes: total #nodes: bias modelled as unit with constant activity)
* *isbias : array[1..numberOfNodes] set 1 if node is bias else 0
* *nodeFirst : array[1..numberOfNodes] first node connected to this unit
* *nodeLast: : array[1..numberOfNodes] last node connected to this unit
* *wFirst : array[1..numberOfNodes] first index in *w for this unit
* *wLast : array[1..numberOfNodes] last (inclusive the bias)
* LEARNING:
* int *wSelected : array[1..nWeights] weights selected for minimization
* double *deriv : array[1..nNodes] derivative of nonlinearity at node
* *error : array[1..nNodes] the error at node
* *dw : array[1..nWeights] total derivative for weights
* *dwi : array[1..nWeights] derivative per pattern
* integer dimension : dimension of minimizer space (<= my nWeights)
* integer nPatterns : the #patterns to be learned
* double **inputPattern: matrix[1..nPatterns][1..nInputs]
* double **targetActivation: matrix[1..nPatterns][1..nOutputs]
* int *wSelected : array[1..numberOfWeights] weights selected for minimization
* double *deriv : array[1..numberOfNodes] derivative of nonlinearity at node
* *error : array[1..numberOfNodes] the error at node
* *dw : array[1..numberOfWeights] total derivative for weights
* *dwi : array[1..numberOfWeights] derivative per pattern
* integer dimension : dimension of minimizer space (<= my numberOfWeights)
* integer numberOfPatterns : the #patterns to be learned
* MAT inputPattern: matrix[1..numberOfPatterns][1..numberOfInputs]
* MAT targetActivation: matrix[1..numberOfPatterns][1..numberOfOutputs]
* double accumulatedCost : accumulated costs of testing/training with patterns
*
* A network consists of nLayers layers. Layer numbering is from 0...nLayers.
* Layer 0 is the input layer, the highest numbered layer is the output layer
* (nLayers <= 4)
* A network consists of numberOfLayers layers. Layer numbering is from 1...numberOfLayers.
* The highest numbered layer is the output layer
* (numberOfLayers <= 4)
* Each layer consists of a number of units. The biases of all the units in a layer
* are modelled with connections to an extra unit in the lower layer (with constant
* activity 1.0). Nodes refers to 'units' + 'bias units'.
* The variable 'nNodes' is the total number of nodes (inclusive bias nodes).
* The variable 'numberOfNodes' is the total number of nodes (inclusive bias nodes).
* E.g. the topology (2,3,4), i.e., 2 inputs, 3 units in the first layer
* and 4 units in the second layer (outputs) is modelled
* with (2+1)+ (3+1)+ (4) = 11 nodes.
* The numbering of the weights is as follows (indices 1..nWeights):
* The numbering of the weights is as follows (indices 1..numberOfWeights):
* E.g., topology (I,H,O) (I inputs, H hidden units and O output units)
* There are a total of H* (I+1) + O* (H+1) weights in this net.
* w[1] - w[I] : I (1)->H (1), I (2)->H (1) ... I (I)->H (1)
......@@ -98,17 +98,17 @@
*
* A number of auxiliary arrays for efficient calculations have been setup.
* For a node k we need to know:
* 1. isbias[1..nNodes] : usage: if (isbias[k]) ...
* true if node k is a bias node. There are nLayers bias nodes
* 2. nodeFirst[1..nNodes] : usage is j=nodeFirst[k];
* 1. isbias[1..numberOfNodes] : usage: if (isbias[k]) ...
* true if node k is a bias node. There are numberOfLayers bias nodes
* 2. nodeFirst[1..numberOfNodes] : usage is j=nodeFirst[k];
* j is the first node that is connected to k .
* 3. nodeLast[1..nNodes] : usage is j=nodeLast[k]
* 3. nodeLast[1..numberOfNodes] : usage is j=nodeLast[k]
* j is the last node that is connected to k (bias included).
* For the calculation of the errors, during learning, in unit k we need to
* know which weights from the preceeding layer connect to it.
* 4. wFirst[1..nNodes] : usage j=wFirst[k]
* 4. wFirst[1..numberOfNodes] : usage j=wFirst[k]
* w[j] is first weight to node k.
* 5. wLast[1..nNodes] : usage j=wLast[k]
* 5. wLast[1..numberOfNodes] : usage j=wLast[k]
* w[j] is last weight to node k.
*/
......@@ -156,16 +156,16 @@ const char32* FFNet_getCategoryOfOutputUnit (FFNet me, integer outputUnit);
integer FFNet_getOutputUnitOfCategory (FFNet me, const char32* category);
void FFNet_propagateToLayer (FFNet me, const double input[], double activity[], integer layer);
void FFNet_propagateToLayer (FFNet me, constVEC input, VEC activity, integer layer);
/* propagate the input through the net to layer and calculate the activities */
void FFNet_propagate (FFNet me, const double input[], double output[]);
void FFNet_propagate (FFNet me, constVEC input, autoVEC *output);
/* step (1) feed forward input from "input layer" to "output layer"
* if output != nullptr the output activity is copied into output.
* postcondition: my activities defined
*/
double FFNet_computeError (FFNet me, const double target[]);
double FFNet_computeError (FFNet me, constVEC target);
/* step (2) calculate error on output nodes w.r.t. desired output */
/* step (3) backpropagate this error to previous nodes */
/* precondition: step (1) */
......@@ -186,9 +186,9 @@ integer FFNet_dimensionOfSearchSpace (FFNet me);
/* count the selected weights */
integer FFNet_getNumberOfWeights (FFNet me);
/* return my nWeights */
/* return my numberOfWeights */
void FFNet_weightConnectsUnits (FFNet me, integer index, integer *fromUnit, integer *toUnit, integer *layer);
void FFNet_weightConnectsUnits (FFNet me, integer index, integer *out_fromUnit, integer *out_toUnit, integer *out_layer);
/*
* w[index] connects unit fromUnit in "layer-1" with unit toUnit in "layer".
* fromUnit returns 0 then w[index] is bias.
......@@ -196,14 +196,11 @@ void FFNet_weightConnectsUnits (FFNet me, integer index, integer *fromUnit, inte
integer FFNet_getNodeNumberFromUnitNumber (FFNet me, integer unit, integer layer);
void FFNet_nodeToUnitInLayer (FFNet me, integer node, integer *unit, integer *layer);
/* translate node index to unit "unit" in layer "layer" */
integer FFNet_getNumberOfLayers (FFNet me);
integer FFNet_getNumberOfUnits (FFNet me);
integer FFNet_getNumberOfHiddenLayers (FFNet me);
integer FFNet_getNumberOfHiddenumberOfLayers (FFNet me);
integer FFNet_getNumberOfUnitsInLayer (FFNet me, int layer);
......
FFNet/FFNet_ActivationList_Categories.cpp
View file @ afb7171d
/* FFNet_ActivationList_Categories.cpp
*
* Copyright (C) 1997-2011, 2015-2017 David Weenink
* Copyright (C) 1997-2011, 2015-2018 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
......@@ -25,43 +25,41 @@
#include "FFNet_ActivationList_Categories.h"
static integer winnerTakesAll (FFNet me, const double activation[]) {
static integer winnerTakesAll (FFNet me, constVEC activation) {
integer pos = 1;
double max = activation[1];
for (integer i = 2; i <= my nOutputs; i ++) {
if (activation [i] > max) {
max = activation [i]; pos = i;
}
for (integer i = 2; i <= my numberOfOutputs; i ++) {
if (activation [i] > max)
max = activation [pos = i];
}
return pos;
}
static integer stochastic (FFNet me, const double activation []) {
static integer stochastic (FFNet me, constVEC activation) {
integer i;
double range = 0.0, lower = 0.0;
for (i = 1; i <= my nOutputs; i ++) {
for (i = 1; i <= my numberOfOutputs; i ++)
range += activation [i];
}
double number = NUMrandomUniform (0.0, range);
for (i = 1; i <= my nOutputs; i ++) {
for (i = 1; i <= my numberOfOutputs; i ++) {
lower += activation [i];
if (number < lower) {
break;
}
if (number < lower) break;
}
return i;
}
autoCategories FFNet_ActivationList_to_Categories (FFNet me, ActivationList activation, int labeling) {
try {
integer (*labelingFunction) (FFNet me, const double act []);
Melder_require (my outputCategories, U"No Categories (has the FFNet been trained yet?).");
Melder_require (my nOutputs == activation -> nx, U"Number of columns and number of outputs should be equal.");
integer (*labelingFunction) (FFNet me, constVEC act);
Melder_require (my outputCategories,
U"No Categories (has the FFNet been trained yet?).");
Melder_require (my numberOfOutputs == activation -> nx,
U"Number of columns and number of outputs should be equal.");
autoCategories thee = Categories_create ();
labelingFunction = labeling == 2 ? stochastic : winnerTakesAll;
for (integer i = 1; i <= activation->ny; i ++) {
integer index = labelingFunction (me, activation -> z [i]);
integer index = labelingFunction (me, activation -> z.row (i));
autoSimpleString item = Data_copy (my outputCategories->at [index]);
thy addItem_move (item.move());
}
......@@ -74,18 +72,19 @@ autoCategories FFNet_ActivationList_to_Categories (FFNet me, ActivationList acti
autoActivationList FFNet_Categories_to_ActivationList (FFNet me, Categories thee) {
try {
autoCategories uniq = Categories_selectUniqueItems (thee);
Melder_require (my outputCategories, U"The FFNet does not have categories.");
Melder_require (my outputCategories,
U"The FFNet does not have categories.");
integer nl = OrderedOfString_isSubsetOf (uniq.get(), my outputCategories.get(), 0);
Melder_require (nl > 0, U"The Categories should match the categories of the FFNet.");
Melder_require (nl > 0,
U"The Categories should match the categories of the FFNet.");
autoActivationList him = ActivationList_create (thy size, my nOutputs);
autoActivationList him = ActivationList_create (thy size, my numberOfOutputs);
for (integer i = 1; i <= thy size; i ++) {
SimpleString category = thy at [i];
integer pos = OrderedOfString_indexOfItem_c (my outputCategories.get(), category -> string.get());
if (pos < 1) {
if (pos < 1)
Melder_throw (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 @ afb7171d
/* FFNet_Eigen.cpp
*
* Copyright (C) 1994-2011, 2015 David Weenink
* Copyright (C) 1994-2018 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
......@@ -30,14 +30,14 @@ void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pc
integer numberOfEigenvalues = eigen -> numberOfEigenvalues;
integer dimension = eigen -> dimension;
if (ix > numberOfEigenvalues || iy > numberOfEigenvalues || my nInputs != dimension) {
if (ix > numberOfEigenvalues || iy > numberOfEigenvalues || my numberOfInputs != dimension)
return;
}
Melder_assert (ix > 0 && iy > 0);
double x1, x2, y1, y2;
if (xmax <= xmin || ymax <= ymin) {
if (xmax <= xmin || ymax <= ymin)
Graphics_inqWindow (g, & x1, & x2, & y1, & y2);
}
if (xmax <= xmin) {
xmin = x1;
xmax = x2;
......@@ -48,12 +48,12 @@ void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pc
}
Graphics_setInner (g);
Graphics_setWindow (g, xmin, xmax, ymin, ymax);
for (integer i = 1; i <= my nUnitsInLayer [1]; i ++) {
integer unitOffset = my nInputs + 1;
for (integer i = 1; i <= my numberOfUnitsInLayer [1]; i ++) {
integer unitOffset = my numberOfInputs + 1;
double c1 = 0.0, c2 = 0.0, bias = my w [my wLast [unitOffset + i]];
double x [6], y [6], xs [3], ys [3];
integer ns = 0;
for (integer j = 1; j <= my nInputs; j ++) {
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];
}
......@@ -63,21 +63,19 @@ void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pc
double p1 = c1 * x [j] + c2 * y [j] + bias;
double p2 = c1 * x [j + 1] + c2 * y [j + 1] + bias;
double r = fabs (p1) / (fabs (p1) + fabs (p2));
if (p1 *p2 > 0 || r == 0.0) {
if (p1 *p2 > 0 || r == 0.0)
continue;
}
if (++ ns > 2) {
break;
}
if (++ ns > 2) break;
xs [ns] = x [j] + (x [j + 1] - x [j]) * r;
ys [ns] = y [j] + (y [j + 1] - y [j]) * r;
}
if (ns < 2) {
if (ns < 2)
Melder_casual (U"Intersection for unit ", i, U" outside range");
} else {
else
Graphics_line (g, xs [1], ys [1], xs [2], ys [2]);
}
}
Graphics_unsetInner (g);
}
......@@ -88,34 +86,28 @@ void FFNet_Eigen_drawIntersection (FFNet me, Eigen eigen, Graphics g, integer pc
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)
{
if (layer < 1 || layer > my nLayers) {
return;
}
if (unit < 1 || unit > my nUnitsInLayer [layer]) {
return;
}
if (pcx > thy numberOfEigenvalues || pcy > thy numberOfEigenvalues) {
return;
}
if (my nUnitsInLayer [layer - 1] != thy dimension) {
return;
}
if (layer < 1 || layer > my numberOfLayers) return;
if (unit < 1 || unit > my numberOfUnitsInLayer [layer]) return;
if (pcx > thy numberOfEigenvalues || pcy > thy numberOfEigenvalues) return;
integer numberOfUnitsInLayer_m1 = ( layer == 1 ? my numberOfInputs : my numberOfUnitsInLayer [layer - 1] );
if (numberOfUnitsInLayer_m1 != thy dimension) return;
double x1, x2, y1, y2;
Graphics_inqWindow (g, & x1, & x2, & y1, & y2);
if (xmax <= xmin) {
xmin = x1; xmax = x2;
xmin = x1;
xmax = x2;
}
if (ymax <= ymin) {
ymin = y1; ymax = y2;
ymin = y1;
ymax = y2;
}
Graphics_setInner (g);
Graphics_setWindow (g, xmin, xmax, ymin, ymax);
integer node = FFNet_getNodeNumberFromUnitNumber (me, unit, layer);
if (node < 1) {
return;
}
if (node < 1) return;
/*
Suppose p1 and p2 are the two points in the eigenplane, spanned by the eigenvectors
......@@ -145,7 +137,7 @@ void FFNet_Eigen_drawDecisionPlaneInEigenspace (FFNet me, Eigen thee, Graphics g
integer iw = my wFirst [node] - 1;
double we1 = 0.0, we2 = 0.0;
for (integer i = 1; i <= my nUnitsInLayer [layer - 1]; i ++) {
for (integer i = 1; i <= numberOfUnitsInLayer_m1; i ++) {
we1 += my w [iw + i] * thy eigenvectors [pcx] [i];
we2 += my w [iw + i] * thy eigenvectors [pcy] [i];
}
......@@ -154,10 +146,12 @@ void FFNet_Eigen_drawDecisionPlaneInEigenspace (FFNet me, Eigen thee, Graphics g
x1 = xmin; x2 = xmax;
y1 = ymin; y2 = ymax;
if (we1 != 0.0) {
x1 = -bias / we1; y1 = 0.0;
x1 = -bias / we1;
y1 = 0.0;
}
if (we2 != 0.0) {
x2 = 0.0; y2 = -bias / we2;
x2 = 0.0;
y2 = -bias / we2;
}
if (we1 == 0.0 && we2 == 0.0) {
Melder_warning (U"We cannot draw the intersection of the neural net decision plane\n"
......@@ -166,11 +160,11 @@ void FFNet_Eigen_drawDecisionPlaneInEigenspace (FFNet me, Eigen thee, Graphics g
}
double xi [3], yi [3]; /* Intersections */
double ni = NUMgetIntersectionsWithRectangle (x1, y1, x2, y2, xmin, ymin, xmax, ymax, xi, yi);
if (ni == 2) {
if (ni == 2)
Graphics_line (g, xi [1], yi [1], xi [2], yi [2]);
} else {
else
Melder_warning (U"There were no intersections in the drawing area.\nPlease enlarge the drawing area.");
}
Graphics_unsetInner (g);
}
......
FFNet/FFNet_Matrix.cpp
View file @ afb7171d
/* FFNet_Matrix.cpp
*
* Copyright (C) 1997-2017 David Weenink
* Copyright (C) 1997-2018 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
......@@ -25,20 +25,20 @@
autoMatrix FFNet_weightsToMatrix (FFNet me, integer layer, bool deltaWeights) {
try {
Melder_require (layer > 0 && layer <= my nLayers, U"Layer should be in [1, ", my nLayers, U"].");
Melder_require (layer > 0 && layer <= my numberOfLayers, U"Layer should be in [1, ", my numberOfLayers, U"].");
integer numberOfUnitsInPreviousLayer = ( layer == 1 ? my numberOfInputs : my numberOfUnitsInLayer [layer - 1] );
autoMatrix thee = Matrix_create (0.5, my nUnitsInLayer [layer] + 0.5, my nUnitsInLayer [layer], 1.0, 1.0,
0.5, my nUnitsInLayer [layer - 1] + 1 + 0.5, my nUnitsInLayer [layer - 1] + 1, 1.0, 1.0);
integer node = 1;
for (integer i = 0; i < layer; i ++) {
node += my nUnitsInLayer [i] + 1;
}
for (integer i = 1; i <= my nUnitsInLayer [layer]; i ++, node ++) {
autoMatrix thee = Matrix_create (0.5, my numberOfUnitsInLayer [layer] + 0.5, my numberOfUnitsInLayer [layer], 1.0, 1.0,
0.5, numberOfUnitsInPreviousLayer + 1 + 0.5, numberOfUnitsInPreviousLayer + 1, 1.0, 1.0);
integer node = my numberOfInputs + 1 + 1;
for (integer i = 1; i < layer; i ++)
node += my numberOfUnitsInLayer [i] + 1;
for (integer i = 1; i <= my numberOfUnitsInLayer [layer]; i ++, node ++) {
integer k = 1;
for (integer j = my wFirst [node]; j <= my wLast [node]; j ++) {
for (integer j = my wFirst [node]; j <= my wLast [node]; j ++)
thy z [k ++] [i] = deltaWeights ? my dwi [j] : my w [j];
}
}
return thee;
} catch (MelderError) {
Melder_throw (me, U": no Matrix created.");
......@@ -47,23 +47,25 @@ autoMatrix FFNet_weightsToMatrix (FFNet me, integer layer, bool deltaWeights) {
autoFFNet FFNet_weightsFromMatrix (FFNet me, Matrix him, integer layer) {
try {
Melder_require (layer > 0 && layer <= my nLayers, U"Layer should be in [1, ", my nLayers, U"].");
Melder_require (my nUnitsInLayer [layer] == his nx, U"The number of columns (", his nx, U") should equal the number of units (", my nUnitsInLayer [layer], U") in layer ", layer, U".");
Melder_require (layer > 0 && layer <= my numberOfLayers,
U"Layer should be in [1, ", my numberOfLayers, U"].");
Melder_require (my numberOfUnitsInLayer [layer] == his nx,
U"The number of columns (", his nx, U") should equal the number of units (", my numberOfUnitsInLayer [layer], U") in layer ", layer, U".");
integer nunits = my nUnitsInLayer [layer - 1] + 1;
Melder_require (nunits == his ny, U"The number of rows (", his ny, U") should equal the number of units (", nunits , U") in layer ", layer - 1, U".");
integer nunits = ( layer == 1 ? my numberOfInputs + 1 : my numberOfUnitsInLayer [layer - 1] + 1 );
Melder_require (nunits == his ny,
U"The number of rows (", his ny, U") should equal the number of units (", nunits , U") in layer ", layer - 1, U".");
autoFFNet thee = Data_copy (me);
integer node = 1;
for (integer i = 0; i < layer; i ++) {
node += thy nUnitsInLayer [i] + 1;
}
for (integer i = 1; i <= thy nUnitsInLayer [layer]; i ++, node ++) {
integer node = my numberOfInputs + 1 + 1;
for (integer i = 1; i < layer; i ++)
node += thy numberOfUnitsInLayer [i] + 1;
for (integer i = 1; i <= thy numberOfUnitsInLayer [layer]; i ++, node ++) {
integer k = 1;
for (integer j = thy wFirst [node]; j <= thy wLast [node]; j ++, k ++) {
for (integer j = thy wFirst [node]; j <= thy wLast [node]; j ++, k ++)
thy w [j] = his z [k] [i];
}
}
return thee;
} catch (MelderError) {
Melder_throw (me, U": no FFNet created.");
......
FFNet/FFNet_Matrix.h
View file @ afb7171d
......@@ -29,8 +29,8 @@
/* The Matrix organization is as follows: */
/* */
/* nx = nUnitsInLayer[layer] */
/* ny = nUnitsInLayer[layer-1]+1 */
/* nx = numberOfUnitsInLayer[layer] */
/* ny = numberOfUnitsInLayer[layer-1]+1 */
/* xmin = 1 xmax = nx */
/* ymin = 1 ymax = ny */
/* dx = dy = 1 */
......
FFNet/FFNet_PatternList.cpp
View file @ afb7171d
/* FFNet_PatternList.cpp
*
* Copyright (C) 1997-2011, 2016 David Weenink
* Copyright (C) 1997-2018 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,10 +24,9 @@
#include "FFNet_PatternList.h"
void FFNet_PatternList_drawActivation (FFNet me, PatternList pattern, Graphics g, integer index) {
if (index < 1 || index > pattern->ny) {
if (index < 1 || index > pattern->ny)
return;
}
FFNet_propagate (me, pattern->z[index], nullptr);
FFNet_propagate (me, pattern -> z.row (index), nullptr);
FFNet_drawActivation (me, g);
}
......
FFNet/FFNet_PatternList_ActivationList.cpp
View file @ afb7171d
/* FFNet_PatternList_ActivationList.cpp
*
* Copyright (C) 1994-2017 David Weenink, 2015,2017 Paul Boersma
* Copyright (C) 1994-2018 David Weenink, 2015,2017 Paul Boersma
*
* 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,52 +27,45 @@
#include "Graphics.h"
#include "FFNet_PatternList_ActivationList.h"
static double func (Daata object, const double p []) {
static double func (Daata object, VEC p) {
FFNet me = (FFNet) object;
Minimizer thee = my minimizer.get();
double fp = 0.0;
longdouble fp = 0.0;
for (integer j = 1, k = 1; k <= my nWeights; k ++) {
for (integer j = 1, k = 1; k <= my numberOfWeights; k ++) {
my dw [k] = 0.0;
if (my wSelected [k]) {
if (my wSelected [k])
my w [k] = p [j ++];
}
}
for (integer i = 1; i <= my nPatterns; i ++) {
FFNet_propagate (me, my inputPattern [i], nullptr);
fp += FFNet_computeError (me, my targetActivation [i]);
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 nWeights; k ++) {
for (integer k = 1; k <= my numberOfWeights; k ++)
my dw [k] += my dwi [k];
}
}
thy funcCalls ++;
return fp;
return (double) fp;
}
static void dfunc_optimized (Daata object, const double p [], double dp []) {
static void dfunc_optimized (Daata object, VEC /* p */, VEC dp) {
FFNet me = (FFNet) object;
(void) p;
integer j = 1;
for (integer k = 1; k <= my nWeights; k ++) {
if (my wSelected [k]) {
for (integer k = 1; k <= my numberOfWeights; k ++) {
if (my wSelected [k])
dp [j ++] = my dw [k];
}
}
}
static void _FFNet_PatternList_ActivationList_checkDimensions (FFNet me, PatternList p, ActivationList a) {
if (my nInputs != p -> nx) {
Melder_throw (U"The PatternList and the FFNet do not match.\nThe number of columns in the PatternList must equal the number of inputs in the FFNet.");
}
if (my nOutputs != a -> nx) {
Melder_throw (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.");
}
if (p -> ny != a -> ny) {
Melder_throw (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.");
}
Melder_require (my numberOfInputs == p -> nx,
U"The PatternList and the FFNet do not match.\nThe number of columns in the PatternList must equal the number of inputs in the FFNet.");
Melder_require (my numberOfOutputs == a -> nx,
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)) {
......@@ -86,33 +79,33 @@ static void _FFNet_PatternList_ActivationList_learn (FFNet me, PatternList patte
// Link the things to be learned
my nPatterns = pattern -> ny;
my inputPattern = pattern -> z;
my targetActivation = activation -> z;
my numberOfPatterns = pattern -> ny;
my inputPattern = pattern -> z.get();
my targetActivation = activation -> z.get();
FFNet_setCostFunction (me, costFunctionType);
if (reset) {
autoNUMvector<double> wbuf (1, my dimension);
autoVEC wbuf = newVECzero (my dimension);
integer k = 1;
for (integer i = 1; i <= my nWeights; i ++) {
for (integer i = 1; i <= my numberOfWeights; i ++) {
if (my wSelected [i]) {
wbuf [k ++] = my w [i];
}
}
Minimizer_reset (my minimizer.get(), wbuf.peek());
Minimizer_reset (my minimizer.get(), wbuf.get());
}
Minimizer_minimize (my minimizer.get(), maxNumOfEpochs, tolerance, 1);
// Unlink
my nPatterns = 0;
my inputPattern = nullptr;
my targetActivation = nullptr;
my numberOfPatterns = 0;
my inputPattern = MAT();
my targetActivation = MAT();
} catch (MelderError) {
my nPatterns = 0;
my inputPattern = nullptr;
my targetActivation = nullptr;
my numberOfPatterns = 0;
my inputPattern = MAT();
my targetActivation = MAT();
}
}
......@@ -158,8 +151,8 @@ double FFNet_PatternList_ActivationList_getCosts_total (FFNet me, PatternList p,
double cost = 0.0;
for (integer i = 1; i <= p -> ny; i ++) {
FFNet_propagate (me, p -> z [i], nullptr);
cost += FFNet_computeError (me, a -> z [i]);
FFNet_propagate (me, p -> z.row (i), nullptr);
cost += FFNet_computeError (me, a -> z.row (i));
}
return cost;
} catch (MelderError) {
......@@ -174,17 +167,19 @@ double FFNet_PatternList_ActivationList_getCosts_average (FFNet me, PatternList
autoActivationList FFNet_PatternList_to_ActivationList (FFNet me, PatternList p, integer layer) {
try {
if (layer < 1 || layer > my nLayers) {
layer = my nLayers;
if (layer < 1 || layer > my numberOfLayers) {
layer = my numberOfLayers;
}
Melder_require (my nInputs == p -> nx, U"The number of colums in the PatternList (", p -> nx, U") should equal the number of inputs in the FFNet (", my nInputs, U").");
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 (my numberOfInputs == p -> nx,
U"The number of colums in the PatternList (", p -> nx, U") should equal the number of inputs in the FFNet (", my numberOfInputs, U").");
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.");
integer nPatterns = p -> ny;
autoActivationList thee = ActivationList_create (nPatterns, my nUnitsInLayer [layer]);
integer numberOfPatterns = p -> ny;
autoActivationList thee = ActivationList_create (numberOfPatterns, my numberOfUnitsInLayer [layer]);
for (integer i = 1; i <= nPatterns; i ++) {
FFNet_propagateToLayer (me, p -> z [i], thy z [i], layer);
for (integer i = 1; i <= numberOfPatterns; i ++) {
FFNet_propagateToLayer (me, p -> z.row (i), thy z.row (i), layer);
}
return thee;
} catch (MelderError) {
......
FFNet/FFNet_PatternList_ActivationList.h
View file @ afb7171d
......@@ -43,6 +43,6 @@ double FFNet_PatternList_ActivationList_getCosts_average (FFNet me, PatternList
autoActivationList FFNet_PatternList_to_ActivationList (FFNet me, PatternList p, integer layer);
/* Calculate the activations at a layer */
/* if (layer<1 || layer > my nLayers) layer = my nLayers; */
/* if (layer<1 || layer > my numberOfLayers) layer = my numberOfLayers; */
#endif /* _FFNet_PatternList_ActivationList_h_ */
FFNet/FFNet_PatternList_Categories.cpp
View file @ afb7171d
/* FFNet_PatternList_Categories.cpp
*
* Copyright (C) 1994-2017 David Weenink
* Copyright (C) 1994-2018 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,9 +28,12 @@
#include "FFNet_PatternList_ActivationList.h"
static void _FFNet_PatternList_Categories_checkDimensions (FFNet me, PatternList p, Categories c) {
Melder_require (my nInputs == p -> nx, U"The PatternList and the FFNet do not match.\nThe number of colums in the PatternList should equal the number of inputs in the FFNet.");
Melder_require (p -> ny == c->size, U"The PatternList and the categories do not match.\nThe number of rows in the PatternList should equal the number of categories.");
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 (my numberOfInputs == p -> nx,
U"The PatternList and the FFNet do not match.\nThe number of colums in the PatternList should equal the number of inputs in the FFNet.");
Melder_require (p -> ny == c->size,
U"The PatternList and the categories do not match.\nThe number of rows in the PatternList should equal the number of categories.");
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.");
}
double FFNet_PatternList_Categories_getCosts_total (FFNet me, PatternList p, Categories c, int costFunctionType) {
......@@ -66,14 +69,17 @@ void FFNet_PatternList_Categories_learnSM (FFNet me, PatternList p, Categories c
autoCategories FFNet_PatternList_to_Categories (FFNet me, PatternList thee, int labeling) {
try {
Melder_require (my outputCategories, U"The FFNet has no output categories.");
Melder_require (my nInputs == thy nx, U"The number of colums in the PatternList (", thy nx, U") should equal the number of inputs in the FFNet (", my nInputs, U").");
Melder_require (_PatternList_checkElements (thee), U"All PatternList elements should be in the interval [0, 1].\nYou could use \"Formula...\" to scale the PatternList values first.");
Melder_require (my outputCategories,
U"The FFNet has no output categories.");
Melder_require (my numberOfInputs == thy nx,
U"The number of colums in the PatternList (", thy nx, U") should equal the number of inputs in the FFNet (", my numberOfInputs, U").");
Melder_require (_PatternList_checkElements (thee),
U"All PatternList elements should be in the interval [0, 1].\nYou could use \"Formula...\" to scale the PatternList values first.");
autoCategories him = Categories_create ();
for (integer k = 1; k <= thy ny; k ++) {
FFNet_propagate (me, thy z [k], nullptr);
FFNet_propagate (me, thy z.row (k), nullptr);
integer index = FFNet_getWinningUnit (me, labeling);
autoSimpleString item = Data_copy (my outputCategories->at [index]);
his addItem_move (item.move());
......
FFNet/FFNet_def.h
View file @ afb7171d
/* FFNet_def.h
*
* Copyright (C) 1994-2008 David Weenink
* Copyright (C) 1994-2018 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
......@@ -26,57 +26,65 @@
#define ooSTRUCT FFNet
oo_DEFINE_CLASS (FFNet, Daata)
oo_INTEGER (nLayers) /* number of layers */
oo_INTEGER_VECTOR_FROM (nUnitsInLayer, 0, nLayers)
oo_INTEGER (numberOfLayers) /* number of layers */
oo_FROM (1)
oo_INTEGER (numberOfInputs)
oo_INTEGER (numberOfOutputs)
oo_ENDFROM
#if oo_READING
oo_VERSION_UNTIL (1)
oo_INTVEC (numberOfUnitsInLayer, numberOfLayers + 1)
numberOfInputs = numberOfUnitsInLayer [1];
numberOfOutputs = numberOfUnitsInLayer [numberOfLayers + 1];
for (integer ilayer = 1; ilayer <= numberOfLayers; ilayer ++)
numberOfUnitsInLayer [ilayer] = numberOfUnitsInLayer [ilayer + 1];
numberOfUnitsInLayer.resize (numberOfLayers);
oo_VERSION_ELSE
oo_INTVEC (numberOfUnitsInLayer, numberOfLayers)
oo_VERSION_END
#else
oo_INTVEC (numberOfUnitsInLayer, numberOfLayers)
#endif
oo_INT (outputsAreLinear)
oo_INT (nonLinearityType)
oo_INT (costFunctionType)
oo_COLLECTION (Categories, outputCategories, SimpleString, 0)
oo_INTEGER (nWeights) /* number of weights */
oo_DOUBLE_VECTOR (w, nWeights)
oo_INTEGER (numberOfWeights) /* number of weights */
oo_VEC (w, numberOfWeights)
#if ! oo_READING && ! oo_WRITING && ! oo_COMPARING
oo_INTEGER (nNodes)
oo_INTEGER (nInputs)
oo_INTEGER (nOutputs)
oo_INTEGER (numberOfNodes)
oo_INTEGER (dimension)
#if oo_DECLARING
double (*nonLinearity) (FFNet /* me */, double /* x */, double * /* deriv */);
void *nlClosure;
double (*costFunction) (FFNet /* me */, const double * /* target */);
double (*costFunction) (FFNet /* me */, constVEC& /* target */);
void *cfClosure;
#endif
#if oo_DECLARING
oo_DOUBLE (accumulatedCost)
oo_INTEGER (nPatterns)
oo_INTEGER (numberOfPatterns)
oo_INTEGER (currentPattern)
double **inputPattern, **targetActivation;
MAT inputPattern, targetActivation;
#endif
#if oo_DECLARING || oo_DESTROYING
oo_OBJECT (Minimizer, 0, minimizer)
#endif
oo_DOUBLE_VECTOR (activity, nNodes)
oo_INTEGER_VECTOR (isbias, nNodes)
oo_INTEGER_VECTOR (nodeFirst, nNodes)
oo_INTEGER_VECTOR (nodeLast, nNodes)
oo_INTEGER_VECTOR (wFirst, nNodes)
oo_INTEGER_VECTOR (wLast, nNodes)
oo_DOUBLE_VECTOR (deriv, nNodes)
oo_DOUBLE_VECTOR (error, nNodes)
oo_INTEGER_VECTOR (wSelected, nWeights)
oo_DOUBLE_VECTOR (dw, nWeights)
oo_DOUBLE_VECTOR (dwi, nWeights)
oo_VEC (activity, numberOfNodes)
oo_INTVEC (isbias, numberOfNodes)
oo_INTVEC (nodeFirst, numberOfNodes)
oo_INTVEC (nodeLast, numberOfNodes)
oo_INTVEC (wFirst, numberOfNodes)
oo_INTVEC (wLast, numberOfNodes)
oo_VEC (deriv, numberOfNodes)
oo_VEC (error, numberOfNodes)
oo_INTVEC (wSelected, numberOfWeights)
oo_VEC (dw, numberOfWeights)
oo_VEC (dwi, numberOfWeights)
#endif
#if oo_READING
......
FFNet/praat_FFNet_init.cpp
View file @ afb7171d
......@@ -168,13 +168,13 @@ DO
DIRECT (INTEGER_FFNet_getNumberOfLayers) {
INTEGER_ONE (FFNet)
integer result = my nLayers;
INTEGER_ONE_END (U" layer", (my nLayers > 1 ? U"s" : U""))
integer result = my numberOfLayers;
INTEGER_ONE_END (U" layer", (my numberOfLayers > 1 ? U"s" : U""))
}
DIRECT (INTEGER_FFNet_getNumberOfOutputs) {
INTEGER_ONE (FFNet)
integer result = my nUnitsInLayer[my nLayers];
integer result = my numberOfUnitsInLayer[my numberOfLayers];
INTEGER_ONE_END (U" units")
}
......@@ -184,13 +184,13 @@ FORM (INTEGER_FFNet_getNumberOfHiddenUnits, U"FFNet: Get number of hidden units"
OK
DO
INTEGER_ONE (FFNet)
integer result = layer > 0 && layer <= my nLayers - 1 ? my nUnitsInLayer[layer] : 0;
integer result = layer > 0 && layer <= my numberOfLayers - 1 ? my numberOfUnitsInLayer[layer] : 0;
INTEGER_ONE_END (U" units")
}
DIRECT (INTEGER_FFNet_getNumberOfInputs) {
INTEGER_ONE (FFNet)
integer result = my nUnitsInLayer[0];
integer result = my numberOfInputs;
INTEGER_ONE_END (U" units")
}
......@@ -199,13 +199,18 @@ FORM (INTEGER_FFNet_getNumberOfHiddenWeights, U"FFNet: Get number of hidden weig
OK
DO
INTEGER_ONE (FFNet)
integer result = (layer > 0 && layer <= my nLayers - 1) ? (my nUnitsInLayer[layer] * (my nUnitsInLayer[layer - 1] + 1)) : 0;
integer result = 0;
if (layer <= my numberOfLayers - 1) {
integer numberOfUnitsInPreviousLayer = ( layer == 1 ? my numberOfInputs : my numberOfUnitsInLayer[layer - 1] );
result = my numberOfUnitsInLayer[layer] * (numberOfUnitsInPreviousLayer + 1);
}
INTEGER_ONE_END (U" weights (including biases)")
}
DIRECT (INTEGER_FFNet_getNumberOfOutputWeights) {
INTEGER_ONE (FFNet)
integer result = my nUnitsInLayer[my nLayers] * (my nUnitsInLayer[my nLayers - 1] + 1);
integer numberOfUnitsInPreviousLayer = ( my numberOfLayers == 1 ? my numberOfInputs : my numberOfUnitsInLayer[my numberOfLayers - 1] );
integer result = my numberOfUnitsInLayer[my numberOfLayers] * (numberOfUnitsInPreviousLayer + 1);
INTEGER_ONE_END (U" weights");
}
......
LPC/Cepstrogram.cpp
View file @ afb7171d
......@@ -40,7 +40,6 @@ autoCepstrogram Cepstrogram_create (double tmin, double tmax, integer nt, double
double qmin, double qmax, integer nq, double dq, double q1) {
try {
autoCepstrogram me = Thing_new (Cepstrogram);
Matrix_init (me.get(), tmin, tmax, nt, dt, t1, qmin, qmax, nq, dq, q1);
return me;
} catch (MelderError) {
......@@ -52,7 +51,6 @@ autoPowerCepstrogram PowerCepstrogram_create (double tmin, double tmax, integer
double qmin, double qmax, integer nq, double dq, double q1) {
try {
autoPowerCepstrogram me = Thing_new (PowerCepstrogram);
Matrix_init (me.get(), tmin, tmax, nt, dt, t1, qmin, qmax, nq, dq, q1);
return me;
} catch (MelderError) {
......@@ -73,14 +71,15 @@ void PowerCepstrogram_paint (PowerCepstrogram me, Graphics g, double tmin, doubl
for (integer i = 1; i <= my ny; i ++) {
for (integer j = 1; j <= my nx; j ++) {
double val = TO10LOG (my z [i] [j]);
min = val < min ? val : min;
max = val > max ? val : max;
if (val < min) min = val;
if (val > max) max = val;
thy z [i] [j] = val;
}
}
double dBminimum = dBmaximum - dynamicRangedB;
if (autoscaling) {
dBminimum = min; dBmaximum = max;
dBminimum = min;
dBmaximum = max;
}
for (integer j = 1; j <= my nx; j ++) {
......@@ -91,18 +90,14 @@ void PowerCepstrogram_paint (PowerCepstrogram me, Graphics g, double tmin, doubl
}
}
double factor = dynamicCompression * (max - lmax);
for (integer i = 1; i <= my ny; i ++) {
thy z [i] [j] += factor;
}
thy z.column(j) += factor;
}
Graphics_setInner (g);
Graphics_setWindow (g, tmin, tmax, qmin, qmax);
Graphics_image (g, thy z,
itmin, itmax,
Graphics_image (g, thy z.part (ifmin, ifmax, itmin, itmax),
Matrix_columnToX (thee.get(), itmin - 0.5),
Matrix_columnToX (thee.get(), itmax + 0.5),
ifmin, ifmax,
Matrix_rowToY (thee.get(), ifmin - 0.5),
Matrix_rowToY (thee.get(), ifmax + 0.5),
dBminimum, dBmaximum);
......@@ -121,13 +116,9 @@ void PowerCepstrogram_subtractTilt_inplace (PowerCepstrogram me, double qstartFi
try {
autoPowerCepstrum thee = PowerCepstrum_create (my ymax, my ny);
for (integer i = 1; i <= my nx; i ++) {
for (integer j = 1; j <= my ny; j ++) {
thy z [1] [j] = my z [j] [i];
}
thy z.row (1) <<= my z.column (i);
PowerCepstrum_subtractTilt_inplace (thee.get(), qstartFit, qendFit, lineType, fitMethod);
for (integer j = 1; j <= my ny; j ++) {
my z [j] [i] = thy z [1] [j];
}
my z.column (i) <<= thy z.row (1);
}
} catch (MelderError) {
Melder_throw (me, U": no tilt subtracted (inline).");
......@@ -150,9 +141,7 @@ autoTable PowerCepstrogram_to_Table_hillenbrand (PowerCepstrogram me, double pit
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 ++) {
for (integer j = 1; j <= my ny; j ++) {
his z [1] [j] = my z [j] [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);
......@@ -171,9 +160,7 @@ autoTable PowerCepstrogram_to_Table_cpp (PowerCepstrogram me, double pitchFloor,
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 ++) {
for (integer j = 1; j <= my ny; j ++) {
his z [1] [j] = my z [j] [i];
}
his z.row (1) <<= my z.column (i);
double qpeak, cpp = PowerCepstrum_getPeakProminence (him.get(), pitchFloor, pitchCeiling, interpolation,
qstartFit, qendFit, lineType, fitMethod, & qpeak);
double rnr = PowerCepstrum_getRNR (him.get(), pitchFloor, pitchCeiling, deltaF0);
......@@ -196,35 +183,23 @@ autoPowerCepstrogram PowerCepstrogram_smooth (PowerCepstrogram me, double timeAv
// 1. average across time
integer numberOfFrames = Melder_ifloor (timeAveragingWindow / my dx);
if (numberOfFrames > 1) {
autoNUMvector<double> qin (1, my nx);
autoNUMvector<double> qout (1, my nx);
autoVEC qin = newVECraw (my nx);
autoVEC qout = newVECraw (my nx);
for (integer iq = 1; iq <= my ny; iq ++) {
for (integer iframe = 1; iframe <= my nx; iframe ++) {
//qin[iframe] = TO10LOG (my z[iq][iframe]);
qin [iframe] = thy z [iq] [iframe];
}
NUMvector_smoothByMovingAverage (qin.peek(), my nx, numberOfFrames, qout.peek());
for (integer iframe = 1; iframe <= my nx; iframe ++) {
//thy z[iq][iframe] = FROMLOG (qout[iframe]); // inverse
thy z [iq] [iframe] = qout [iframe]; // inverse
}
qin.all() <<= thy z.row (iq); // ppgb: why this extra copying?
VECsmoothByMovingAverage_preallocated (qout.get(), qin.get(), numberOfFrames);
thy z.row (iq) <<= qout.all();
}
}
// 2. average across quefrencies
integer numberOfQuefrencyBins = Melder_ifloor (quefrencyAveragingWindow / my dy);
if (numberOfQuefrencyBins > 1) {
autoNUMvector<double> qin (1, thy ny);
autoNUMvector<double> qout (1, thy ny);
autoVEC qin = newVECraw (thy ny);
autoVEC qout = newVECraw (thy ny);
for (integer iframe = 1; iframe <= my nx; iframe ++) {
for (integer iq = 1; iq <= thy ny; iq ++) {
//qin[iq] = TO10LOG (my z[iq][iframe]);
qin [iq] = thy z [iq] [iframe];
}
NUMvector_smoothByMovingAverage (qin.peek(), thy ny, numberOfQuefrencyBins, qout.peek());
for (integer iq = 1; iq <= thy ny; iq ++) {
//thy z[iq][iframe] = FROMLOG (qout[iq]);
thy z [iq] [iframe] = qout [iq];
}
qin.get() <<= thy z.column (iframe);
VECsmoothByMovingAverage_preallocated (qout.get(), qin.get(), numberOfQuefrencyBins);
thy z.column (iframe) <<= qout.get();
}
}
return thee;
......@@ -248,9 +223,7 @@ autoPowerCepstrum PowerCepstrogram_to_PowerCepstrum_slice (PowerCepstrogram me,
integer iframe = Sampled_xToNearestIndex (me, time);
iframe = iframe < 1 ? 1 : iframe > my nx ? my nx : iframe;
autoPowerCepstrum thee = PowerCepstrum_create (my ymax, my ny);
for (integer i = 1; i <= my ny; i ++) {
thy z [1] [i] = my z [i] [iframe];
}
thy z.row (1) <<= my z.column (iframe);
return thee;
} catch (MelderError) {
Melder_throw (me, U": Cepstrum not extracted.");
......@@ -300,14 +273,13 @@ autoPowerCepstrogram Sound_to_PowerCepstrogram (Sound me, double pitchFloor, dou
Sounds_multiply (sframe.get(), window.get());
autoSpectrum spec = Sound_to_Spectrum (sframe.get(), true); // FFT yes
autoPowerCepstrum cepstrum = Spectrum_to_PowerCepstrum (spec.get());
for (integer i = 1; i <= nq; i ++) {
for (integer i = 1; i <= nq; i ++)
thy z [i] [iframe] = cepstrum -> z [1] [i];
}
if ((iframe % 10) == 1) {
if ((iframe % 10) == 1)
Melder_progress ((double) iframe / nFrames, U"PowerCepstrogram analysis of frame ",
iframe, U" out of ", nFrames, U".");
}
}
return thee;
} catch (MelderError) {
Melder_throw (me, U": no PowerCepstrogram created.");
......@@ -319,18 +291,16 @@ autoCepstrum Spectrum_to_Cepstrum_hillenbrand (Spectrum me) {
try {
autoNUMfft_Table fftTable;
// originalNumberOfSamplesProbablyOdd irrelevant
if (my x1 != 0.0) {
Melder_throw (U"A Fourier-transformable Spectrum must have a first frequency of 0 Hz, not ", my x1, U" Hz.");
}
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;
autoCepstrum thee = Cepstrum_create (0.5 / my dx, my nx);
NUMfft_Table_init (& fftTable, my nx);
autoNUMvector<double> amp (1, my nx);
autoVEC amp = newVECraw (my nx);
for (integer i = 1; i <= my nx; i ++) {
for (integer i = 1; i <= my nx; i ++)
amp [i] = my v_getValueAtSample (i, 0, 2);
}
NUMfft_forward (& fftTable, amp.peek());
NUMfft_forward (& fftTable, amp.get());
for (integer i = 1; i <= my nx; i ++) {
double val = amp [i] / numberOfSamples;// scaling 1/n because ifft(fft(1))= n;
......@@ -346,18 +316,18 @@ autoCepstrum Spectrum_to_Cepstrum_hillenbrand (Spectrum me) {
// re im re im re im
// ((fft[1],0) (fft[2],fft[3]), (,), (,), (fft[nfft], 0)) nfft even
// ((fft[1],0) (fft[2],fft[3]), (,), (,), (fft[nfft-1], fft[nfft])) nfft uneven
static void complexfftoutput_to_power (double *fft, integer nfft, double *dbs, bool to_db) {
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 = (nfft + 2) / 2;
integer nend = nfft % 2 == 0 ? nfftdiv2p1 : nfftdiv2p1 + 1;
integer nfftdiv2p1 = (fft.size + 2) / 2;
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];
valsq = re * re + im * im;
dbs [i] = to_db ? TOLOG (valsq) : valsq;
}
if (nfft % 2 == 0) {
valsq = fft [nfft] * fft [nfft];
if (fft.size % 2 == 0) {
valsq = fft [fft.size] * fft [fft.size];
dbs[nfftdiv2p1] = to_db ? TOLOG (valsq) : valsq;
}
}
......@@ -365,69 +335,54 @@ static void complexfftoutput_to_power (double *fft, integer nfft, double *dbs, b
autoPowerCepstrogram Sound_to_PowerCepstrogram_hillenbrand (Sound me, double pitchFloor, double dt) {
try {
// minimum analysis window has 3 periods of lowest pitch
double analysisWidth = 3.0 / pitchFloor;
if (analysisWidth > my dx * my nx) {
analysisWidth = my dx * my nx;
}
const double physicalDuration = my dx * my nx;
const double analysisWidth = std::min (3.0 / pitchFloor, physicalDuration);
double samplingFrequency = 1.0 / my dx;
autoSound thee;
if (samplingFrequency > 30000) {
if (samplingFrequency > 30000.0) {
samplingFrequency = samplingFrequency / 2.0;
thee = Sound_resample (me, samplingFrequency, 1);
} else {
thee = Data_copy (me);
}
// pre-emphasis with fixed coefficient 0.9
for (integer i = thy nx; i > 1; i --) {
for (integer i = thy nx; i > 1; i --)
thy z [1] [i] -= 0.9 * thy z [1] [i - 1];
}
integer nosInWindow = (integer) floorl (analysisWidth * samplingFrequency), numberOfFrames;
if (nosInWindow < 8) {
Melder_throw (U"Analysis window too short.");
}
integer nosInWindow = Melder_ifloor (analysisWidth * samplingFrequency), numberOfFrames;
Melder_require (nosInWindow >= 8, U"Analysis window too short.");
double t1;
Sampled_shortTermAnalysis (thee.get(), analysisWidth, dt, & numberOfFrames, & t1);
autoNUMvector<double> hamming (1, nosInWindow);
for (integer i = 1; i <= nosInWindow; i ++) {
hamming [i] = 0.54 -0.46 * cos(2 * NUMpi * (i - 1) / (nosInWindow - 1));
}
autoVEC hamming = newVECraw (nosInWindow);
for (integer i = 1; i <= nosInWindow; i ++)
hamming [i] = 0.54 - 0.46 * cos (2.0 * NUMpi * (i - 1) / (nosInWindow - 1));
integer nfft = 8; // minimum possible
while (nfft < nosInWindow) { nfft *= 2; }
integer nfftdiv2 = nfft / 2;
autoNUMvector<double> fftbuf (1, nfft); // "complex" array
autoNUMvector<double> spectrum (1, nfftdiv2 + 1); // +1 needed
const integer nfftdiv2 = nfft / 2;
autoVEC fftbuf = newVECzero (nfft); // "complex" array
autoVEC spectrum = newVECzero (nfftdiv2 + 1); // +1 needed
autoNUMfft_Table fftTable;
NUMfft_Table_init (& fftTable, nfft); // sound to spectrum
double qmax = 0.5 * nfft / samplingFrequency, dq = qmax / (nfftdiv2 + 1);
const double qmax = 0.5 * nfft / samplingFrequency, dq = qmax / (nfftdiv2 + 1);
autoPowerCepstrogram him = PowerCepstrogram_create (my xmin, my xmax, numberOfFrames, dt, t1, 0, qmax, nfftdiv2+1, dq, 0);
autoMelderProgress progress (U"Cepstrogram analysis");
for (integer iframe = 1; iframe <= numberOfFrames; iframe ++) {
double tbegin = t1 + (iframe - 1) * dt - analysisWidth / 2;
tbegin = tbegin < thy xmin ? thy xmin : tbegin;
integer istart = Sampled_xToLowIndex (thee.get(), tbegin); // ppgb: afronding naar beneden?
istart = istart < 1 ? 1 : istart;
integer iend = istart + nosInWindow - 1;
iend = iend > thy nx ? thy nx : iend;
for (integer i = 1; i <= nosInWindow; i ++) {
fftbuf [i] = thy z [1] [istart + i - 1] * hamming [i];
}
for (integer i = nosInWindow + 1; i <= nfft; i ++) {
fftbuf [i] = 0;
}
NUMfft_forward (&fftTable, fftbuf.peek());
complexfftoutput_to_power (fftbuf.peek(), nfft, spectrum.peek(), true); // log10(|fft|^2)
// subtract average
double specmean = spectrum [1];
for (integer i = 2; i <= nfftdiv2 + 1; i ++) {
specmean += spectrum [i];
}
specmean /= nfftdiv2 + 1;
for (integer i = 1; i <= nfftdiv2 + 1; i ++) {
spectrum [i] -= specmean;
}
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
fftbuf.part (1, nosInWindow) <<= thy z.row (1).part (istart, iend) * hamming.all();
fftbuf.part (nosInWindow + 1, nfft) <<= 0.0;
NUMfft_forward (& fftTable, fftbuf.get());
complexfftoutput_to_power (fftbuf.get(), spectrum.get(), true); // log10(|fft|^2)
VECcentre_inplace (spectrum.get()); // subtract average
/*
* Here we diverge from Hillenbrand as he takes the fft of half of the spectral values.
* H. forgets that the actual spectrum has nfft/2+1 values. Thefore, we take the inverse
......@@ -438,18 +393,17 @@ autoPowerCepstrogram Sound_to_PowerCepstrogram_hillenbrand (Sound me, double pit
fftbuf [1] = spectrum [1];
for (integer i = 2; i < nfftdiv2 + 1; i ++) {
fftbuf [i+i-2] = spectrum [i];
fftbuf [i+i-1] = 0;
fftbuf [i+i-1] = 0.0;
}
fftbuf [nfft] = spectrum [nfftdiv2 + 1];
NUMfft_backward (& fftTable, fftbuf.peek());
for (integer i = 1; i <= nfftdiv2 + 1; i ++) {
NUMfft_backward (& fftTable, fftbuf.get());
for (integer i = 1; i <= nfftdiv2 + 1; i ++)
his z [i] [iframe] = fftbuf [i] * fftbuf [i];
}
if ((iframe % 10) == 1) {
if (iframe % 10 == 1)
Melder_progress ((double) iframe / numberOfFrames, U"Cepstrogram analysis of frame ",
iframe, U" out of ", numberOfFrames, U".");
}
}
return him;
} catch (MelderError) {
Melder_throw (me, U": no Cepstrogram created.");
......@@ -459,9 +413,9 @@ autoPowerCepstrogram Sound_to_PowerCepstrogram_hillenbrand (Sound me, double pit
double PowerCepstrogram_getCPPS (PowerCepstrogram me, bool subtractTiltBeforeSmoothing, double timeAveragingWindow, double quefrencyAveragingWindow, double pitchFloor, double pitchCeiling, double deltaF0, int interpolation, double qstartFit, double qendFit, int lineType, int fitMethod) {
try {
autoPowerCepstrogram flattened;
if (subtractTiltBeforeSmoothing) {
if (subtractTiltBeforeSmoothing)
flattened = PowerCepstrogram_subtractTilt (me, qstartFit, qendFit, lineType, fitMethod);
}
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);
......@@ -474,9 +428,9 @@ double PowerCepstrogram_getCPPS (PowerCepstrogram me, bool subtractTiltBeforeSmo
double PowerCepstrogram_getCPPS_hillenbrand (PowerCepstrogram me, bool subtractTiltBeforeSmoothing, double timeAveragingWindow, double quefrencyAveragingWindow, double pitchFloor, double pitchCeiling) {
try {
autoPowerCepstrogram him;
if (subtractTiltBeforeSmoothing) {
if (subtractTiltBeforeSmoothing)
him = PowerCepstrogram_subtractTilt (me, 0.001, 0, 1, 1);
}
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);
......
LPC/Cepstrum.cpp
View file @ afb7171d
......@@ -92,32 +92,28 @@ static void _Cepstrum_draw (Cepstrum me, Graphics g, double qmin, double qmax, d
Graphics_setInner (g);
if (qmax <= qmin) {
qmin = my xmin; qmax = my xmax;
qmin = my xmin;
qmax = my xmax;
}
integer imin, imax;
if (! Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax)) {
return;
}
autoNUMvector<double> y (imin, imax);
integer numberOfSelected = imax - imin + 1;
autoVEC y = newVECraw (numberOfSelected);
for (integer i = imin; i <= imax; i ++) {
y [i] = my v_getValueAtSample (i, (power ? 1 : 0), 0);
for (integer i = 1; i <= numberOfSelected; i ++) {
y [i] = my v_getValueAtSample (imin + i - 1, (power ? 1 : 0), 0);
}
if (autoscaling) {
NUMvector_extrema (y.peek(), imin, imax, & minimum, & maximum);
} else {
for (integer i = imin; i <= imax; i ++) {
if (y [i] > maximum) {
y [i] = maximum;
} else if (y [i] < minimum) {
y [i] = minimum;
}
}
}
if (autoscaling)
NUMextrema (y.get(), & minimum, & maximum);
else
VECclip_inplace (y.get(), minimum, maximum);
Graphics_setWindow (g, qmin, qmax, minimum, maximum);
Graphics_function (g, y.peek(), imin, imax, Matrix_columnToX (me, imin), Matrix_columnToX (me, imax));
Graphics_function (g, y.at, 1, numberOfSelected, Matrix_columnToX (me, imin), Matrix_columnToX (me, imax));
Graphics_unsetInner (g);
......@@ -143,7 +139,8 @@ void PowerCepstrum_drawTiltLine (PowerCepstrum me, Graphics g, double qmin, doub
Graphics_setInner (g);
if (qmax <= qmin) {
qmin = my xmin; qmax = my xmax;
qmin = my xmin;
qmax = my xmax;
}
if (dBminimum >= dBmaximum) { // autoscaling
......@@ -164,7 +161,8 @@ void PowerCepstrum_drawTiltLine (PowerCepstrum me, Graphics g, double qmin, doub
Graphics_setWindow (g, qmin, qmax, dBminimum, dBmaximum);
qend = qend == 0 ? my xmax : qend;
if (qend <= qstart) {
qend = my xmax; qstart = my xmin;
qend = my xmax;
qstart = my xmin;
}
qstart = qstart < my xmin ? my xmin : qstart;
qend = qend > my xmax ? my xmax : qend;
......@@ -184,13 +182,13 @@ void PowerCepstrum_drawTiltLine (PowerCepstrum me, Graphics g, double qmin, doub
qstart = 0.1 * dq; // some small offset to avoid log(0)
n--;
}
autoNUMvector<double> y (1, n);
autoVEC y = newVECraw (n);
for (integer i = 1; i <= n; i ++) {
double q = q1 + (i - 1) * dq;
y [i] = a * log (q) + intercept;
}
Graphics_function (g, y.peek(), 1, n, qstart, qend);
Graphics_function (g, y.at, 1, n, qstart, qend);
} else {
double y1 = a * qstart + intercept, y2 = a * qend + intercept;
if (y1 >= dBminimum && y2 >= dBminimum) {
......@@ -213,35 +211,33 @@ void PowerCepstrum_drawTiltLine (PowerCepstrum me, Graphics g, double qmin, doub
* method == 1 : Least squares fit
* method == 2 : Theil's partial robust fit
*/
void PowerCepstrum_fitTiltLine (PowerCepstrum me, double qmin, double qmax, double *p_a, double *p_intercept, int lineType, int method) {
void PowerCepstrum_fitTiltLine (PowerCepstrum me, double qmin, double qmax, double *out_a, double *out_intercept, int lineType, int method) {
try {
double a, intercept;
if (qmax <= qmin) {
qmin = my xmin; qmax = my xmax;
qmin = my xmin;
qmax = my xmax;
}
integer imin, imax;
if (! Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax)) {
if (! Matrix_getWindowSamplesX (me, qmin, qmax, & imin, & imax))
return;
}
imin = (lineType == 2 && imin == 1) ? 2 : imin; // log(0) is undefined!
integer numberOfPoints = imax - imin + 1;
if (numberOfPoints < 2) {
Melder_throw (U"Not enough points for fit.");
}
autoNUMvector<double> y (1, numberOfPoints);
autoNUMvector<double> x (1, numberOfPoints);
Melder_require (numberOfPoints > 1, U"Not enough points for fit.");
autoVEC y = newVECraw (numberOfPoints);
autoVEC x = newVECraw (numberOfPoints);
for (integer i = 1; i <= numberOfPoints; i ++) {
integer isamp = imin + i - 1;
x [i] = my x1 + (isamp - 1) * my dx;
if (lineType == 2) {
if (lineType == 2)
x [i] = log (x [i]);
}
y [i] = my v_getValueAtSample (isamp, 1, 0);
}
if (method == 3) { // try local maxima first
autoNUMvector<double> ym (1, numberOfPoints / 2 + 1);
autoNUMvector<double> xm (1, numberOfPoints / 2 + 1);
autoVEC ym = newVECraw (numberOfPoints / 2 + 1);
autoVEC xm = newVECraw (numberOfPoints / 2 + 1);
integer numberOfLocalPeaks = 0;
// forget y [1] if y [2]<y [1] and y [n] if y [n-1]<y [n] !
for (integer i = 2; i <= numberOfPoints; i ++) {
......@@ -252,16 +248,17 @@ void PowerCepstrum_fitTiltLine (PowerCepstrum me, double qmin, double qmax, doub
}
if (numberOfLocalPeaks > numberOfPoints / 10) {
for (integer i = 1; i <= numberOfLocalPeaks; i ++) {
x [i] = xm [i]; y [i] = ym [i];
x [i] = xm [i];
y [i] = ym [i];
}
numberOfPoints = numberOfLocalPeaks;
}
method = 2; // robust fit of peaks
}
// fit a straight line through (x,y)'s
NUMlineFit (x.peek(), y.peek(), numberOfPoints, & a, & intercept, method);
if (p_intercept) { *p_intercept = intercept; }
if (p_a) { *p_a = a; }
NUMlineFit (x.get(), y.get(), & a, & intercept, method);
if (out_intercept) *out_intercept = intercept;
if (out_a) *out_a = a;
} catch (MelderError) {
Melder_throw (me, U": couldn't fit a line.");
}
......@@ -319,20 +316,14 @@ void PowerCepstrum_smooth_inplace (PowerCepstrum me, double quefrencyAveragingWi
try {
integer numberOfQuefrencyBins = Melder_ifloor (quefrencyAveragingWindow / my dx);
if (numberOfQuefrencyBins > 1) {
autoNUMvector<double> qin (1, my nx);
autoNUMvector<double> qout (1, my nx);
for (integer iq = 1; iq <= my nx; iq ++) {
qin [iq] = my z [1] [iq];
}
double *xin, *xout;
for (integer k = 1; k <= numberOfIterations; k ++) {
xin = k % 2 == 1 ? qin.peek() : qout.peek ();
xout = k % 2 == 1 ? qout.peek () : qin.peek();
NUMvector_smoothByMovingAverage (xin, my nx, numberOfQuefrencyBins, xout);
}
for (integer iq = 1; iq <= my nx; iq ++) {
my z [1] [iq] = xout [iq];
}
autoVEC qin = newVECcopy (my z.row (1));
autoVEC qout = newVECraw (my nx);
for (integer k = 1; k <= numberOfIterations; k ++)
if (k % 2 == 1)
VECsmoothByMovingAverage_preallocated (qout.get(), qin.get(), numberOfQuefrencyBins);
else
VECsmoothByMovingAverage_preallocated (qin.get(), qout.get(), numberOfQuefrencyBins);
my z.row (1) <<= ( numberOfIterations % 2 == 1 ? qout.get() : qin.get() );
}
} catch (MelderError) {
Melder_throw (me, U": not smoothed.");
......@@ -345,62 +336,37 @@ autoPowerCepstrum PowerCepstrum_smooth (PowerCepstrum me, double quefrencyAverag
return thee;
}
void PowerCepstrum_getMaximumAndQuefrency (PowerCepstrum me, double pitchFloor, double pitchCeiling, int interpolation, double *p_peakdB, double *p_quefrency) {
double peakdB, quefrency;
void PowerCepstrum_getMaximumAndQuefrency (PowerCepstrum me, double pitchFloor, double pitchCeiling, int interpolation, double *out_peakdB, double *out_quefrency) {
autoPowerCepstrum thee = Data_copy (me);
double lowestQuefrency = 1.0 / pitchCeiling, highestQuefrency = 1.0 / pitchFloor;
for (integer i = 1; i <= my nx; i ++) {
thy z [1] [i] = my v_getValueAtSample (i, 1, 0); // 10 log val^2
}
double peakdB, quefrency;
Vector_getMaximumAndX ((Vector) thee.get(), lowestQuefrency, highestQuefrency, 1, interpolation, & peakdB, & quefrency); // FIXME cast
if (p_peakdB) {
*p_peakdB = peakdB;
}
if (p_quefrency) {
*p_quefrency = quefrency;
}
}
#if 0
static void Cepstrum_getZ (Cepstrum me, integer imin, integer imax, double peakdB, double slope, double intercept, int lineType, double *z) {
integer ndata = imax - imin + 1;
autoNUMvector<double> dabs (1, ndata);
for (integer i = imin; i <= imax; i ++) {
double q = my x1 + (i - 1) * my dx;
q = ( i == 1 ? 0.5 * my dx : q ); // approximation
double xq = ( lineType == 2 ? log (q) : q );
double db_background = slope * xq + intercept;
double db_cepstrum = my v_getValueAtSample (i, 1, 0);
double diff = exp ((db_cepstrum - db_background) * NUMln10 / 10.0) - 1e-30;
//double diff = fabs (db_cepstrum - db_background);
dabs [i - imin + 1] = diff;
}
double q50 = NUMquantile (ndata, dabs.peek(), 0.5);
double peak = exp (peakdB * NUMln10 / 10.0) - 1e-30;
if (z) {
*z = peak / q50;
if (out_peakdB) *out_peakdB = peakdB;
if (out_quefrency) *out_quefrency = quefrency;
}
}
#endif
double PowerCepstrum_getRNR (PowerCepstrum me, double pitchFloor, double pitchCeiling, double f0fractionalWidth) {
double rnr = undefined;
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))
return rnr;
}
integer ndata = imax - imin + 1;
if (ndata < 2) {
if (ndata < 2)
return rnr;
}
// how many peaks in interval ?
integer npeaks = 2;
while (qpeak > 0 && qpeak * npeaks <= qmax) { npeaks++; }
while (qpeak > 0 && qpeak * npeaks <= qmax)
npeaks++;
npeaks--;
double sum = 0, sumr = 0;
double sum = 0.0, sumr = 0.0;
for (integer i = imin; i <= imax; i ++) {
double val = my v_getValueAtSample (i, 0, 0);
double qx = my x1 + (i - 1) * my dx;
......@@ -418,32 +384,28 @@ double PowerCepstrum_getRNR (PowerCepstrum me, double pitchFloor, double pitchCe
}
}
}
rnr = sumr >= sum ? 1000000 : sumr / (sum - sumr);
rnr = sumr >= sum ? 1000000.0 : sumr / (sum - sumr);
return rnr;
}
double PowerCepstrum_getPeakProminence_hillenbrand (PowerCepstrum me, double pitchFloor, double pitchCeiling, double *qpeak) {
double PowerCepstrum_getPeakProminence_hillenbrand (PowerCepstrum me, double pitchFloor, double pitchCeiling, double *out_qpeak) {
double slope, intercept, quefrency, peakdB;
PowerCepstrum_fitTiltLine (me, 0.001, 0, & slope, & intercept, 1, 1);
autoPowerCepstrum thee = Data_copy (me);
PowerCepstrum_subtractTiltLine_inplace (thee.get(), slope, intercept, 1);
PowerCepstrum_getMaximumAndQuefrency (thee.get(), pitchFloor, pitchCeiling, 0, & peakdB, & quefrency);
if (qpeak) {
*qpeak = quefrency;
}
if (out_qpeak) *out_qpeak = quefrency;
return peakdB;
}
double PowerCepstrum_getPeakProminence (PowerCepstrum me, double pitchFloor, double pitchCeiling, int interpolation, double qstartFit, double qendFit, int lineType, int fitMethod, double *p_qpeak) {
double PowerCepstrum_getPeakProminence (PowerCepstrum me, double pitchFloor, double pitchCeiling, int interpolation, double qstartFit, double qendFit, int lineType, int fitMethod, double *out_qpeak) {
double slope, intercept, qpeak, peakdB;
PowerCepstrum_fitTiltLine (me, qstartFit, qendFit, &slope, &intercept, lineType, fitMethod);
PowerCepstrum_getMaximumAndQuefrency (me, pitchFloor, pitchCeiling, interpolation, & peakdB, & qpeak);
double xq = lineType == 2 ? log(qpeak) : qpeak;
double db_background = slope * xq + intercept;
double cpp = peakdB - db_background;
if (p_qpeak) {
*p_qpeak = qpeak;
}
if (out_qpeak) *out_qpeak = qpeak;
return cpp;
}
......@@ -473,10 +435,9 @@ autoPowerCepstrum Matrix_to_PowerCepstrum (Matrix me) {
autoPowerCepstrum Matrix_to_PowerCepstrum_row (Matrix me, integer row) {
try {
autoPowerCepstrum thee = PowerCepstrum_create (my xmax, my nx);
if (row < 1 || row > my ny) {
Melder_throw (U"Row number should be between 1 and ", my ny, U" inclusive.");
}
NUMvector_copyElements (my z [row], thy z [1], 1, my nx);
Melder_require (row > 0 && row <= my ny,
U"Row number should be between 1 and ", my ny, U" inclusive.");
thy z.row (1) <<= my z.row (row);
return thee;
} catch (MelderError) {
Melder_throw (me, U": no PowerCepstrum created.");
......@@ -486,12 +447,9 @@ autoPowerCepstrum Matrix_to_PowerCepstrum_row (Matrix me, integer row) {
autoPowerCepstrum Matrix_to_PowerCepstrum_column (Matrix me, integer col) {
try {
autoPowerCepstrum thee = PowerCepstrum_create (my ymax, my ny);
if (col < 1 || col > my nx) {
Melder_throw (U"Column number should be between 1 and ", my nx, U" inclusive.");
}
for (integer i = 1; i <= my ny; i ++) {
thy z [1] [i] = my z [i] [col];
}
Melder_require (col > 0 && col <= my nx,
U"Column number should be between 1 and ", my nx, U" inclusive.");
thy z.row (1) <<= my z.column (col);
return thee;
} catch (MelderError) {
Melder_throw (me, U": no PowerCepstrum created.");
......
LPC/Cepstrum.h
View file @ afb7171d
......@@ -2,7 +2,7 @@
#define _Cepstrum_h_
/* Cepstrum.h
*
* Copyright (C) 1994-2013, 2015-2016 David Weenink
* Copyright (C) 1994-2018 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
......@@ -95,7 +95,7 @@ double PowerCepstrum_getPeakProminence_hillenbrand (PowerCepstrum me, double pit
double PowerCepstrum_getRNR (PowerCepstrum me, double pitchFloor, double pitchCeiling, double f0fractionalWidth);
double PowerCepstrum_getPeakProminence (PowerCepstrum me, double pitchFloor, double pitchCeiling, int interpolation, double qstartFit, double qendFit, int lineType, int fitMethod, double *qpeak);
void PowerCepstrum_fitTiltLine (PowerCepstrum me, double qmin, double qmax, double *slope, double *intercept, int lineType, int method);
void PowerCepstrum_fitTiltLine (PowerCepstrum me, double qmin, double qmax, double *out_slope, double *out_intercept, int lineType, int method);
autoPowerCepstrum PowerCepstrum_subtractTilt (PowerCepstrum me, double qstartFit, double qendFit, int lineType, int fitMethod);
void PowerCepstrum_subtractTilt_inplace (PowerCepstrum me, double qstartFit, double qendFit, int lineType, int fitMethod);
......