diff --git a/find_songs_from_hand_annotations.m b/find_songs_from_hand_annotations.m new file mode 100644 index 0000000..89f1c92 --- /dev/null +++ b/find_songs_from_hand_annotations.m @@ -0,0 +1,100 @@ +function likelihoodModels = find_songs_from_hand_annotations(male_songs,female_songs,overlap_songs) +%Inputs arrays of selected male, female, and overlap song and outputs +%likelihood models +%Inputs: +% male_songs -> N_male x length(frequencies) array of male song wavelet amplitudes +% female_songs -> N_female x length(frequencies) array of female song wavelet amplitudes +% overlap_songs -> N_overlap x length(frequencies) array of overlapping song amplitudes +% +%Output: +% likelihoodModels -> struct containing likelihood models needed to run +% segmentVirilisSong.m +% +% (C) Gordon J. Berman, Jan Clemens, Kelly M. LaRue, and Mala Murthy, 2015 +% Princeton University + + + %load parameters + segParams = params_virilis(); + + probModes = segParams.probModes; + maxNumPeaks = segParams.maxNumPeaks; + maxNumPeaks_firstMode = segParams.maxNumPeaks_firstMode; + frequencies = segParams.fc; + gmm_replicates = segParams.gmm_replicates; + maxNumGMM = segParams.maxNumGMM; + + fprintf(1,'Finding Male Principal Components\n'); + male_mean = mean(male_songs); + [coeffs_male,scores_male,latent_male] = princomp(male_songs); + + + fprintf(1,'Finding Female Principal Components\n'); + female_mean = mean(female_songs); + [coeffs_female,scores_female,latent_female] = princomp(female_songs); + + + if nargin >= 3 && isempty(overlap_songs) + fprintf(1,'Finding Overlap Principal Components\n'); + both_mean = mean(overlap_songs); + [coeffs_both,scores_both,latent_both] = princomp(overlap_songs); + end + + + + fprintf(1,'Finding PDFs\n'); + + malePDFs = cell(probModes,1); + femalePDFs = cell(probModes,1); + if nargin >= 3 && isempty(overlap_songs) + bothPDFs = cell(probModes,1); + end + for i=1:probModes + + fprintf(1,'\t #%2i of %2i\n',i,probModes); + if i == 1 + q = maxNumPeaks_firstMode; + else + q = maxNumPeaks; + end + + malePDFs{i} = findBestGMM_AIC(scores_male(:,i),q,gmm_replicates,maxNumGMM); + femalePDFs{i} = findBestGMM_AIC(scores_female(:,i),q,gmm_replicates,maxNumGMM); + if nargin >= 3 && isempty(overlap_songs) + bothPDFs{i} = findBestGMM_AIC(scores_both(:,i),q,gmm_replicates,maxNumGMM); + end + end + + + likelihoodModels.malePDFs = malePDFs; + likelihoodModels.femalePDFs = femalePDFs; + likelihoodModels.noisePDFs = noisePDFs; + + likelihoodModels.male_mean = male_mean; + likelihoodModels.female_mean = female_mean; + likelihoodModels.noise_mean = noise_mean; + + likelihoodModels.coeffs_male = coeffs_male; + likelihoodModels.coeffs_female = coeffs_female; + likelihoodModels.coeffs_noise = coeffs_noise; + + likelihoodModels.latent_male = latent_male; + likelihoodModels.latent_female = latent_female; + likelihoodModels.latent_noise = latent_noise; + + likelihoodModels.probModes = probModes; + likelihoodModels.frequencies = frequencies; + + K = scal2frq(1,wav,segParams.dt); + scales = K ./ frequencies; + likelihoodModels.scales = scales; + + + if nargin >= 3 && isempty(overlap_songs) + likelihoodModels.bothPDFs = bothPDFs; + likelihoodModels.both_mean = both_mean; + likelihoodModels.coeffs_both = coeffs_both; + likelihoodModels.latent_both = latent_both; + end + + diff --git a/find_songs_from_hand_annotations.m~ b/find_songs_from_hand_annotations.m~ new file mode 100644 index 0000000..a876b33 --- /dev/null +++ b/find_songs_from_hand_annotations.m~ @@ -0,0 +1,84 @@ +function likelihoodModels = find_songs_from_hand_annotations(male_songs,female_songs,overlap_songs) +%Inputs arrays of selected male, female, and overlap song, + + %load parameters + segParams = params_virilis(); + + probModes = segParams.probModes; + maxNumPeaks = segParams.maxNumPeaks; + maxNumPeaks_firstMode = segParams.maxNumPeaks_firstMode; + frequencies = segParams.fc; + gmm_replicates = segParams.gmm_replicates; + maxNumGMM = segParams.maxNumGMM; + + fprintf(1,'Finding Male Principal Components\n'); + male_mean = mean(male_songs); + [coeffs_male,scores_male,latent_male] = princomp(male_songs); + + + fprintf(1,'Finding Female Principal Components\n'); + female_mean = mean(female_songs); + [coeffs_female,scores_female,latent_female] = princomp(female_songs); + + + if nargin >= 3 && isempty(overlap_songs) + fprintf(1,'Finding Overlap Principal Components\n'); + both_mean = mean(overlap_songs); + [coeffs_both,scores_both,latent_both] = princomp(overlap_songs); + end + + + + fprintf(1,'Finding PDFs\n'); + + malePDFs = cell(probModes,1); + femalePDFs = cell(probModes,1); + if nargin >= 3 && isempty(overlap_songs) + bothPDFs = cell(probModes,1); + end + for i=1:probModes + + fprintf(1,'\t #%2i of %2i\n',i,probModes); + if i == 1 + q = maxNumPeaks_firstMode; + else + q = maxNumPeaks; + end + + malePDFs{i} = findBestGMM_AIC(scores_male(:,i),q,gmm_replicates,maxNumGMM); + femalePDFs{i} = findBestGMM_AIC(scores_female(:,i),q,gmm_replicates,maxNumGMM); + if nargin >= 3 && isempty(overlap_songs) + bothPDFs{i} = findBestGMM_AIC(scores_both(:,i),q,gmm_replicates,maxNumGMM); + end + end + + + likelihoodModels.malePDFs = malePDFs; + likelihoodModels.femalePDFs = femalePDFs; + likelihoodModels.noisePDFs = noisePDFs; + + likelihoodModels.male_mean = male_mean; + likelihoodModels.female_mean = female_mean; + likelihoodModels.noise_mean = noise_mean; + + likelihoodModels.coeffs_male = coeffs_male; + likelihoodModels.coeffs_female = coeffs_female; + likelihoodModels.coeffs_noise = coeffs_noise; + + likelihoodModels.latent_male = latent_male; + likelihoodModels.latent_female = latent_female; + likelihoodModels.latent_noise = latent_noise; + + likelihoodModels.probModes = probModes; + likelihoodModels.frequencies = frequencies; + likelihoodModels.scales = scales; + + + if nargin >= 3 && isempty(overlap_songs) + likelihoodModels.bothPDFs = bothPDFs; + likelihoodModels.both_mean = both_mean; + likelihoodModels.coeffs_both = coeffs_both; + likelihoodModels.latent_both = latent_both; + end + + diff --git a/params_virilis.m b/params_virilis.m index c761095..b64d0c7 100644 --- a/params_virilis.m +++ b/params_virilis.m @@ -1,159 +1,178 @@ function segParams = params_virilis(fs) - -if nargin < 1 || isempty(fs) - fs = 10000; -end - -%%%%%%%%%%%%%%ALL USER DEFINED PARAMETERS ARE SET HERE%%%%%%%%%%%%%%%%%%%%% - - -%frequencies to use (don't change unless you want to re-run the likelihood -%models) -fc = 100:20:900; - -%factor for computing window around pulse peak -%(this determines how much of the signal before and after the peak -%is included in the pulse, and sets the parameters w0 and w1.) -pulseWindow = round(fs/25); - -%factor times the mean of xempty - only pulses larger than -%this amplitude are counted as true pulses -noiseFactor = 1; - -%male pulse carrier frequency -mpf = 350; - -%100Hz lowpass filter for male pulse detection -male = 100; - -%40Hz lowpass filter for female pulse detection -female = 40; - -%20Hz lowpass filter for male song bout detection -bout = 20; - -%male bout detection based on IPI, in samples -male_IPI = 250; - -%Smoothing parameter for data -filterWindow = 5; - - -%%%%%%%%%%%%%%%%%%%%%%% - -%Threshold for P(both singing | data) -probThreshold = 1; - -%Male probability threshold -maleThreshold = .2; +%Creates parameters data structure +%Input: +% fs -> Data sampling frequency (default = 10000) +% +%Output: +% segParams -> struct containing parameters used in model building and fitting +% +% (C) Gordon J. Berman, Jan Clemens, Kelly M. LaRue, and Mala Murthy, 2015 +% Princeton University + + + if nargin < 1 || isempty(fs) + fs = 10000; + end + + %%%%%%%%%%%%%%ALL USER DEFINED PARAMETERS ARE SET HERE%%%%%%%%%%%%%%%%%%%%% + + + %frequencies to use (don't change unless you want to re-run the likelihood + %models) + fc = 100:20:900; + + %factor for computing window around pulse peak + %(this determines how much of the signal before and after the peak + %is included in the pulse, and sets the parameters w0 and w1.) + pulseWindow = round(fs/25); + + %factor times the mean of xempty - only pulses larger than + %this amplitude are counted as true pulses + noiseFactor = 1; + + %male pulse carrier frequency + mpf = 350; + + %100Hz lowpass filter for male pulse detection + male = 100; + + %40Hz lowpass filter for female pulse detection + female = 40; + + %20Hz lowpass filter for male song bout detection + bout = 20; + + %male bout detection based on IPI, in samples + male_IPI = 250; + + %Smoothing parameter for data + filterWindow = 5; + + + %%%%%%%%%%%%%%%%%%%%%%% + + %Threshold for P(both singing | data) + probThreshold = 1; + + %Male probability threshold + maleThreshold = .2; + + %minimum size of a female pulse in the midst of a male pulse (in data + %points) + minFemalePulseSize = 150; + + %number of PCA modes to use in analysis (3 -> 41) + probModes = 20; %(DONT CHANGE THIS!!!!!!!!!!!!!!) + + %Pnoise < noiseThreshold to count as signal + noiseThreshold = .5; + + %minimum percentage of time where p(male) > p(female) in order for the bout + %to be counted as male (***) + minMaleBoutFraction = 1/5; + + %Number of time points over which to test male activity (i is called male if + %the number of initially called male frames is greater than + %minMaleBoutFraction*maleTestDuration points over + %(i-maleTestDuration):(i+maleTestDuration) + maleTestDuration = 500; + + %minimum duration for a bout to be called a male bout (in time points) (***) + minMaleDuration = 1000; + + %Smoothing parameter for male and combination song likelihoods (***) + smoothParameter_male = 25; + + %Smoothing parameter for female and nosie song likelihoods (***) + smoothParameter_female = 50; + + %Smoothing parameter for female and nosie song likelihoods + smoothParameter_amplitudes = 25; + + %Threshold for amplitude + amplitudeThreshold = .1; + + %Threshold for for noise likelihood (***) + %(smaller = more stringent, do not make lower than -5490, must be changed if probModes is changed) + noiseLikelihoodThreshold = 0; + + %Threshold for calling the posterior value noise (between 0 and 1). + %Closer to 1 implies a more stringent cut-off + ampPostThreshold = .95; + + %Maximum noise data set from which to create noise models + maxNumNoise = 500000; + + %Maximum number of peaks to create GMM noise distribution PDF + maxNumPeaks = 4; + + %Maximum number of peaks to create GMM noise distribution PDF for the first + %mode + maxNumPeaks_firstMode = 6; + + %Number of PDF modes used to find female bouts within male bouts + num_male_both_modes = 5; + + %Definition of closeness for female peak elimination (in ms) + female_IPI_limit = 26; + + %Num of consecutive close female peaks needed for elimination (should be + %odd) + num_female_IPI_limit = 3; + + %Length of time (in seconds) + stop_recording_time = 140; + + %Whether or not to refine probabilities using the PDF projection theorem + usePDFProjections = false; + + %Maximum number of data to use in GMM fitting of projection distributions + maxNumGMM = 10000; + + %Number of replicates to use in GMM fits + gmm_replicates = 3; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + + segParams.maxNumGMM = maxNumGMM; + segParams.gmm_replicates = gmm_replicates; + segParams.fc = fc; + segParams.fs = fs; + segParams.dt = 1/fs; + segParams.pulsewindow = pulseWindow; + segParams.noiseFactor = noiseFactor; + segParams.mpf = mpf; + segParams.male = male; + segParams.female = female; + segParams.bout = bout; + segParams.male_IPI = male_IPI; + segParams.probThreshold = probThreshold; + segParams.minFemalePulseSize = minFemalePulseSize; + segParams.filterWindow = filterWindow; + segParams.probModes = probModes; + segParams.noiseThreshold = noiseThreshold; + segParams.minMaleBoutFraction = minMaleBoutFraction; + segParams.minMaleDuration = minMaleDuration; + segParams.smoothParameter_male = smoothParameter_male; + segParams.smoothParameter_female = smoothParameter_female; + segParams.smoothParameter_amplitudes = smoothParameter_amplitudes; + segParams.amplitudeThreshold = amplitudeThreshold; + segParams.noiseLikelihoodThreshold = noiseLikelihoodThreshold; + segParams.maleThreshold = maleThreshold; + segParams.ampPostThreshold = ampPostThreshold; + segParams.maxNumNoise = maxNumNoise; + segParams.maxNumPeaks = maxNumPeaks; + segParams.maxNumPeaks_firstMode = maxNumPeaks_firstMode; + segParams.num_male_both_modes = num_male_both_modes; + segParams.female_IPI_limit = female_IPI_limit; + segParams.num_female_IPI_limit = num_female_IPI_limit; + segParams.stop_recording_time = stop_recording_time; + segParams.usePDFProjections = usePDFProjections; + segParams.maleTestDuration = maleTestDuration; -%minimum size of a female pulse in the midst of a male pulse (in data -%points) -minFemalePulseSize = 150; - -%number of PCA modes to use in analysis (3 -> 41) -probModes = 20; %(DONT CHANGE THIS!!!!!!!!!!!!!!) - -%Pnoise < noiseThreshold to count as signal -noiseThreshold = .5; - -%minimum percentage of time where p(male) > p(female) in order for the bout -%to be counted as male (***) -minMaleBoutFraction = 1/5; - -%Number of time points over which to test male activity (i is called male if -%the number of initially called male frames is greater than -%minMaleBoutFraction*maleTestDuration points over -%(i-maleTestDuration):(i+maleTestDuration) -maleTestDuration = 500; - -%minimum duration for a bout to be called a male bout (in time points) (***) -minMaleDuration = 1000; - -%Smoothing parameter for male and combination song likelihoods (***) -smoothParameter_male = 25; - -%Smoothing parameter for female and nosie song likelihoods (***) -smoothParameter_female = 50; - -%Smoothing parameter for female and nosie song likelihoods -smoothParameter_amplitudes = 25; - -%Threshold for amplitude -amplitudeThreshold = .1; - -%Threshold for for noise likelihood (***) -%(smaller = more stringent, do not make lower than -5490, must be changed if probModes is changed) -noiseLikelihoodThreshold = 0; - -%Threshold for calling the posterior value noise (between 0 and 1). -%Closer to 1 implies a more stringent cut-off -ampPostThreshold = .95; - -%Maximum noise data set from which to create noise models -maxNumNoise = 500000; - -%Maximum number of peaks to create GMM noise distribution PDF -maxNumPeaks = 4; - -%Maximum number of peaks to create GMM noise distribution PDF for the first -%mode -maxNumPeaks_firstMode = 6; - -%Number of PDF modes used to find female bouts within male bouts -num_male_both_modes = 5; - -%definition of closeness for female peak elimination (in ms) -female_IPI_limit = 26; - -%num of consecutive close female peaks needed for elimination (should be -%odd) -num_female_IPI_limit = 3; - -%length of time (in seconds) -stop_recording_time = 140; - -%whether or not to refine probabilities using the PDF projection theorem -usePDFProjections = false; - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - - - -segParams.fc = fc; -segParams.fs = fs; -segParams.pulsewindow = pulseWindow; -segParams.noiseFactor = noiseFactor; -segParams.mpf = mpf; -segParams.male = male; -segParams.female = female; -segParams.bout = bout; -segParams.male_IPI = male_IPI; -segParams.probThreshold = probThreshold; -segParams.minFemalePulseSize = minFemalePulseSize; -segParams.filterWindow = filterWindow; -segParams.probModes = probModes; -segParams.noiseThreshold = noiseThreshold; -segParams.minMaleBoutFraction = minMaleBoutFraction; -segParams.minMaleDuration = minMaleDuration; -segParams.smoothParameter_male = smoothParameter_male; -segParams.smoothParameter_female = smoothParameter_female; -segParams.smoothParameter_amplitudes = smoothParameter_amplitudes; -segParams.amplitudeThreshold = amplitudeThreshold; -segParams.noiseLikelihoodThreshold = noiseLikelihoodThreshold; -segParams.maleThreshold = maleThreshold; -segParams.ampPostThreshold = ampPostThreshold; -segParams.maxNumNoise = maxNumNoise; -segParams.maxNumPeaks = maxNumPeaks; -segParams.maxNumPeaks_firstMode = maxNumPeaks_firstMode; -segParams.num_male_both_modes = num_male_both_modes; -segParams.female_IPI_limit = female_IPI_limit; -segParams.num_female_IPI_limit = num_female_IPI_limit; -segParams.stop_recording_time = stop_recording_time; -segParams.usePDFProjections = usePDFProjections; -segParams.maleTestDuration = maleTestDuration; - diff --git a/params_virilis.m~ b/params_virilis.m~ new file mode 100644 index 0000000..dfc5398 --- /dev/null +++ b/params_virilis.m~ @@ -0,0 +1,181 @@ +function segParams = params_virilis(fs) +%Fits song data to previously-defined likelihood models +%Inputs: +% data -> 1d array containing a time series of song +% likelihoodModels -> fitted likelihood models from +% find_songs_from_hand_annotation.m (default:exampleLikelihoodModels) +% samplingFrequency -> data sampling frequency in Hz (default = 1e4) +% +%Outputs: +% maleBoutInfo -> struct containing statistis from found male bouts +% femaleBoutInfo -> struct containing statistis from found female bouts +% run_data -> struct containing statistis from analysis +% +% (C) Gordon J. Berman, Jan Clemens, Kelly M. LaRue, and Mala Murthy, 2015 +% Princeton University + + + if nargin < 1 || isempty(fs) + fs = 10000; + end + + %%%%%%%%%%%%%%ALL USER DEFINED PARAMETERS ARE SET HERE%%%%%%%%%%%%%%%%%%%%% + + + %frequencies to use (don't change unless you want to re-run the likelihood + %models) + fc = 100:20:900; + + %factor for computing window around pulse peak + %(this determines how much of the signal before and after the peak + %is included in the pulse, and sets the parameters w0 and w1.) + pulseWindow = round(fs/25); + + %factor times the mean of xempty - only pulses larger than + %this amplitude are counted as true pulses + noiseFactor = 1; + + %male pulse carrier frequency + mpf = 350; + + %100Hz lowpass filter for male pulse detection + male = 100; + + %40Hz lowpass filter for female pulse detection + female = 40; + + %20Hz lowpass filter for male song bout detection + bout = 20; + + %male bout detection based on IPI, in samples + male_IPI = 250; + + %Smoothing parameter for data + filterWindow = 5; + + + %%%%%%%%%%%%%%%%%%%%%%% + + %Threshold for P(both singing | data) + probThreshold = 1; + + %Male probability threshold + maleThreshold = .2; + + %minimum size of a female pulse in the midst of a male pulse (in data + %points) + minFemalePulseSize = 150; + + %number of PCA modes to use in analysis (3 -> 41) + probModes = 20; %(DONT CHANGE THIS!!!!!!!!!!!!!!) + + %Pnoise < noiseThreshold to count as signal + noiseThreshold = .5; + + %minimum percentage of time where p(male) > p(female) in order for the bout + %to be counted as male (***) + minMaleBoutFraction = 1/5; + + %Number of time points over which to test male activity (i is called male if + %the number of initially called male frames is greater than + %minMaleBoutFraction*maleTestDuration points over + %(i-maleTestDuration):(i+maleTestDuration) + maleTestDuration = 500; + + %minimum duration for a bout to be called a male bout (in time points) (***) + minMaleDuration = 1000; + + %Smoothing parameter for male and combination song likelihoods (***) + smoothParameter_male = 25; + + %Smoothing parameter for female and nosie song likelihoods (***) + smoothParameter_female = 50; + + %Smoothing parameter for female and nosie song likelihoods + smoothParameter_amplitudes = 25; + + %Threshold for amplitude + amplitudeThreshold = .1; + + %Threshold for for noise likelihood (***) + %(smaller = more stringent, do not make lower than -5490, must be changed if probModes is changed) + noiseLikelihoodThreshold = 0; + + %Threshold for calling the posterior value noise (between 0 and 1). + %Closer to 1 implies a more stringent cut-off + ampPostThreshold = .95; + + %Maximum noise data set from which to create noise models + maxNumNoise = 500000; + + %Maximum number of peaks to create GMM noise distribution PDF + maxNumPeaks = 4; + + %Maximum number of peaks to create GMM noise distribution PDF for the first + %mode + maxNumPeaks_firstMode = 6; + + %Number of PDF modes used to find female bouts within male bouts + num_male_both_modes = 5; + + %Definition of closeness for female peak elimination (in ms) + female_IPI_limit = 26; + + %Num of consecutive close female peaks needed for elimination (should be + %odd) + num_female_IPI_limit = 3; + + %Length of time (in seconds) + stop_recording_time = 140; + + %Whether or not to refine probabilities using the PDF projection theorem + usePDFProjections = false; + + %Maximum number of data to use in GMM fitting of projection distributions + maxNumGMM = 10000; + + %Number of replicates to use in GMM fits + gmm_replicates = 3; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + + segParams.maxN + segParams.fc = fc; + segParams.fs = fs; + segParams.pulsewindow = pulseWindow; + segParams.noiseFactor = noiseFactor; + segParams.mpf = mpf; + segParams.male = male; + segParams.female = female; + segParams.bout = bout; + segParams.male_IPI = male_IPI; + segParams.probThreshold = probThreshold; + segParams.minFemalePulseSize = minFemalePulseSize; + segParams.filterWindow = filterWindow; + segParams.probModes = probModes; + segParams.noiseThreshold = noiseThreshold; + segParams.minMaleBoutFraction = minMaleBoutFraction; + segParams.minMaleDuration = minMaleDuration; + segParams.smoothParameter_male = smoothParameter_male; + segParams.smoothParameter_female = smoothParameter_female; + segParams.smoothParameter_amplitudes = smoothParameter_amplitudes; + segParams.amplitudeThreshold = amplitudeThreshold; + segParams.noiseLikelihoodThreshold = noiseLikelihoodThreshold; + segParams.maleThreshold = maleThreshold; + segParams.ampPostThreshold = ampPostThreshold; + segParams.maxNumNoise = maxNumNoise; + segParams.maxNumPeaks = maxNumPeaks; + segParams.maxNumPeaks_firstMode = maxNumPeaks_firstMode; + segParams.num_male_both_modes = num_male_both_modes; + segParams.female_IPI_limit = female_IPI_limit; + segParams.num_female_IPI_limit = num_female_IPI_limit; + segParams.stop_recording_time = stop_recording_time; + segParams.usePDFProjections = usePDFProjections; + segParams.maleTestDuration = maleTestDuration; + + diff --git a/segmentVirilisSong.m b/segmentVirilisSong.m index afa5ff6..d358fdb 100644 --- a/segmentVirilisSong.m +++ b/segmentVirilisSong.m @@ -1,4 +1,19 @@ -function [maleBoutInfo,femaleBoutInfo,run_data] = segmentVirilisSong(data,likelihoodModels,samplingFrequency) +function [maleBoutInfo,femaleBoutInfo,run_data] = ... + segmentVirilisSong(data,likelihoodModels,samplingFrequency) +%Fits song data to previously-defined likelihood models +%Inputs: +% data -> 1d array containing a time series of song +% likelihoodModels -> fitted likelihood models from +% find_songs_from_hand_annotation.m (default:exampleLikelihoodModels) +% samplingFrequency -> data sampling frequency in Hz (default = 1e4) +% +%Outputs: +% maleBoutInfo -> struct containing statistis from found male bouts +% femaleBoutInfo -> struct containing statistis from found female bouts +% run_data -> struct containing statistis from analysis +% +% (C) Gordon J. Berman, Jan Clemens, Kelly M. LaRue, and Mala Murthy, 2015 +% Princeton University addpath('utilities'); addpath('subroutines'); diff --git a/subroutines/findProbabilities_wavelet.m b/subroutines/findProbabilities_wavelet.m index 3a263c4..3435f90 100644 --- a/subroutines/findProbabilities_wavelet.m +++ b/subroutines/findProbabilities_wavelet.m @@ -1,7 +1,27 @@ function [probs,likelihoods,noiseP,P,CCs] = ... findProbabilities_wavelet(data,likelihoodModels,noiseModel,... segmentParameters,plotsOn,plotData) - +%Find posterior probabilities for wavelet data from likelihood models +%Inputs: +% data -> 1d time series or num_time_points x length(frequencies) array of wavelet data +% (wavelets performed in the case of the former, but not the latter) +% likelihoodModels -> fitted likelihood models from +% find_songs_from_hand_annotation.m (default:exampleLikelihoodModels) +% noiseModel -> fitted noise likelihood model +% segmentParameters -> struct containing parameters +% plotsOn -> true if plot are desired (default = true) +% plotData -> 1d song data to plot (only used if plotsOn = true) +% +%Outputs: +% probs -> Nx4 array of posterior probabilities +% likelihoods -> Nx4 array of log likelihood scores +% noiseP -> probabilities that a particular point belongs to the noise +% model +% P -> wavelet amplitudes +% CCs -> connected components from male and female data (only if plotsOn = true) +% +% (C) Gordon J. Berman, Jan Clemens, Kelly M. LaRue, and Mala Murthy, 2015 +% Princeton University if nargin < 4 || isempty(plotsOn) plotsOn = true; diff --git a/subroutines/findProbabilities_wavelet.m~ b/subroutines/findProbabilities_wavelet.m~ new file mode 100644 index 0000000..3979d67 --- /dev/null +++ b/subroutines/findProbabilities_wavelet.m~ @@ -0,0 +1,169 @@ +function [probs,likelihoods,noiseP,P,CCs] = ... + findProbabilities_wavelet(data,likelihoodModels,noiseModel,... + segmentParameters,plotsOn,plotData) +%Find posterior probabilities for wavelet data from likelihood models +%Inputs: +% data -> num_time_points x length(frequencies) array of wavelet data +% likelihoodModels -> fitted likelihood models from +% find_songs_from_hand_annotation.m (default:exampleLikelihoodModels) +% noiseModel -> fitted noise likelihood model +% segmentParameters -> struct containing parameters +% plotsOn -> true if plot are desired (default = true) +% plotData -> 1d song data to plot (only used if plotsOn = true) +% +%Outputs: +% probs -> Nx4 array of posterior probabilities +% likelihoods -> Nx4 array of log likelihood scores +% noiseP -> probabilities that a particular point belongs to the noise +% model +% P -> wavelet amplitudes +% CCs -> struct containing statistis from analysis +% +% (C) Gordon J. Berman, Jan Clemens, Kelly M. LaRue, and Mala Murthy, 2015 +% Princeton University + + if nargin < 4 || isempty(plotsOn) + plotsOn = true; + end + + + probModes = segmentParameters.probModes; + + smoothParameter_male = segmentParameters.smoothParameter_male; + smoothParameter_female = segmentParameters.smoothParameter_female; + + min_male_length = 1; + min_female_length = 1; + + malePDFs = likelihoodModels.malePDFs; + femalePDFs = likelihoodModels.femalePDFs; + noisePDFs = noiseModel.noisePDFs; + bothPDFs = likelihoodModels.bothPDFs; + + male_mean = likelihoodModels.male_mean; + female_mean = likelihoodModels.female_mean; + noise_mean = noiseModel.noise_mean; + both_mean = likelihoodModels.both_mean; + + coeffs_male = likelihoodModels.coeffs_male; + coeffs_female = likelihoodModels.coeffs_female; + coeffs_noise = noiseModel.coeffs_noise; + coeffs_both = likelihoodModels.coeffs_both; + + + + s = size(data); + if min(s) == 1 + fprintf(1,'Performing Wavelet Transform\n'); + scales = likelihoodModels.scales; + wav = 'fbsp2-1-2'; + C = cwt(data,scales,wav); + P = C'.*conj(C'); + clear C; + else + P = data; + if plotsOn && nargin == 5 + data = plotData; + end + end + + N = length(P(:,1)); + + + fprintf(1,'Calculating Projections\n'); + dataScores_male = bsxfun(@minus,P,male_mean) * coeffs_male(:,1:probModes); + dataScores_female = bsxfun(@minus,P,female_mean) * coeffs_female(:,1:probModes); + dataScores_noise = bsxfun(@minus,P,noise_mean) * coeffs_noise(:,1:probModes); + dataScores_both = bsxfun(@minus,P,both_mean) * coeffs_both(:,1:probModes); + + + fprintf(1,'Finding Likelihoods\n'); + likelihoods = zeros(N,4); + + for i=1:probModes + + likelihoods(:,1) = likelihoods(:,1) + log(pdf(malePDFs{i},dataScores_male(:,i))); + likelihoods(:,2) = likelihoods(:,2) + log(pdf(femalePDFs{i},dataScores_female(:,i))); + likelihoods(:,3) = likelihoods(:,3) + log(pdf(noisePDFs{i},dataScores_noise(:,i))); + likelihoods(:,4) = likelihoods(:,4) + log(pdf(bothPDFs{i},dataScores_both(:,i))); + + end + + + fprintf(1,'Computing Probabilities\n'); + maxVal = max(likelihoods(~isnan(likelihoods) & ~isinf(likelihoods))); + minVal = min(likelihoods(~isnan(likelihoods) & ~isinf(likelihoods))); + likelihoods(isnan(likelihoods) | isinf(likelihoods)) = minVal; + + if smoothParameter_male > 1 + likelihoods(:,1) = gaussianfilterdata(likelihoods(:,1),smoothParameter_male); + likelihoods(:,4) = gaussianfilterdata(likelihoods(:,4),smoothParameter_male); + end + + if smoothParameter_female > 1 + likelihoods(:,2) = gaussianfilterdata(likelihoods(:,2),smoothParameter_female); + likelihoods(:,3) = gaussianfilterdata(likelihoods(:,3),smoothParameter_female); + end + + probs = exp(likelihoods-maxVal); + partition = sum(probs,2); + probs = bsxfun(@rdivide,probs,partition); + + sumProbs = sum(probs)./N; + + maxVals = max(likelihoods,[],2); + subLikes = exp(bsxfun(@minus,likelihoods,maxVals)); + partition = zeros(N,1); + for i=1:4 + partition = partition + subLikes(:,i)*sumProbs(i); + end + + probs = bsxfun(@rdivide,bsxfun(@times,subLikes,sumProbs),partition); + noiseP = probs(:,3); + + + probs = bsxfun(@rdivide,probs,sum(probs,2)); + + + if plotsOn || nargout == 5 + + probs2 = probs(:,1:3); + probs2(:,1) = probs2(:,1) + probs(:,4); + + [~,maxIdx] = max(probs2,[],2); + + + CC_male = largeBWConnComp(maxIdx == 1 | maxIdx == 4,min_male_length); + CC_female = largeBWConnComp(maxIdx == 2,min_female_length); + + CCs = {CC_male,CC_female}; + + figure + + subplot(2,1,1) + hold on + for i=1:length(CC_male.PixelIdxList) + rectangle('Position',... + [CC_male.PixelIdxList{i}(1) -1 length(CC_male.PixelIdxList{i}) 2],... + 'facecolor','b','edgecolor','b'); + end + + for i=1:length(CC_female.PixelIdxList) + rectangle('Position',... + [CC_female.PixelIdxList{i}(1) -1 length(CC_female.PixelIdxList{i}) 2],... + 'facecolor','r','edgecolor','r'); + end + + plot(data,'k-') + + ylim([-.3 .3]) + + + + subplot(2,1,2) + plot(probs(:,1)+probs(:,4),'bo-'); + hold on + plot(probs(:,2),'rs-'); + plot(probs(:,3),'k^-'); + ylim([-.02 1.02]) + end