    # function [PCG_Features, featuresFs] = getSpringerPCGFeatures(audio_data, Fs, figures)
    
    # Get the features used in the Springer segmentation algorithm. These 
# features include:
# -The homomorphic envelope (as performed in Schmidt et al's paper)
# -The Hilbert envelope
# -A wavelet-based feature
# -A PSD-based feature
# This function was developed for use in the paper:
# D. Springer et al., "Logistic Regression-HSMM-based Heart Sound 
# Segmentation," IEEE Trans. Biomed. Eng., In Press, 2015.
    
    ## INPUTS:
# audio_data: array of data from which to extract features
# Fs: the sampling frequency of the audio data
# figures (optional): boolean variable dictating the display of figures
    
    ## OUTPUTS:
# PCG_Features: array of derived features
# featuresFs: the sampling frequency of the derived features. This is set
# in default_Springer_HSMM_options.py
    
    ## Copyright (C) 2016  David Springer
# dave.springer@gmail.com
# 
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
    
    
@function
def getSpringerPCGFeatures(audio_data=None,Fs=None,figures=None,*args,**kwargs):
    varargin = getSpringerPCGFeatures.varargin
    nargin = getSpringerPCGFeatures.nargin

    # function PCG_Features = getSpringerPCGFeatures(audio, Fs)
# Get the features used in the Springer segmentation algorithm.
    
    if (nargin < 3):
        figures=copy(false)

    
    springer_options=copy(default_Springer_HSMM_options)

    # Check to see if the Wavelet toolbox is available on the machine:
    include_wavelet=springer_options.include_wavelet_feature

    featuresFs=springer_options.audio_segmentation_Fs

    
    ## 25-400Hz 4th order Butterworth band pass
    audio_data=butterworth_low_pass_filter(audio_data,2,400,Fs,false)

    audio_data=butterworth_high_pass_filter(audio_data,2,25,Fs)

    ## Spike removal from the original paper:
    audio_data=schmidt_spike_removal(audio_data,Fs)

    ## Find the homomorphic envelope
    homomorphic_envelope=Homomorphic_Envelope_with_Hilbert(audio_data,Fs)

    # Downsample the envelope:
    downsampled_homomorphic_envelope=resample(homomorphic_envelope,featuresFs,Fs)

    # normalise the envelope:
    downsampled_homomorphic_envelope=normalise_signal(downsampled_homomorphic_envelope)

    ## Hilbert Envelope
    hilbert_envelope=Hilbert_Envelope(audio_data,Fs)

    downsampled_hilbert_envelope=resample(hilbert_envelope,featuresFs,Fs)

    downsampled_hilbert_envelope=normalise_signal(downsampled_hilbert_envelope)

    ## Power spectral density feature:
    
    psd=get_PSD_feature_Springer_HMM(audio_data,Fs,40,60).T

    psd=resample(psd,length(downsampled_homomorphic_envelope),length(psd))

    psd=normalise_signal(psd)

    ## Wavelet features:
    
    if (include_wavelet):
        wavelet_level=3

        #     wavelet_name ='rbio3.9';
        wavelet_name='db7'

        if (length(audio_data) < dot(Fs,1.025)):
            audio_data=concat([[audio_data],[zeros(round(dot(0.025,Fs)),1)]])

        cD,cA=getDWT(audio_data,wavelet_level,wavelet_name,nargout=2)

        wavelet_feature=abs(cD(wavelet_level,arange()))

        wavelet_feature=wavelet_feature(arange(1,length(homomorphic_envelope)))

        downsampled_wavelet=resample(wavelet_feature,featuresFs,Fs)

        downsampled_wavelet=normalise_signal(downsampled_wavelet).T

    
    ##
    
    if (include_wavelet):
        PCG_Features=concat([downsampled_homomorphic_envelope,downsampled_hilbert_envelope,psd,downsampled_wavelet])

    else:
        PCG_Features=concat([downsampled_homomorphic_envelope,downsampled_hilbert_envelope,psd])

    
    ## Plotting figures
    if (figures):
        figure('Name','PCG features')
        t1=(arange(1,length(audio_data))) / Fs

        plot(t1,audio_data)
        hold('on')
        t2=(arange(1,length(PCG_Features))) / featuresFs

        plot(t2,PCG_Features)
        pause()
    