MatLab Functions

MatLab functions for Auditory/Cochlear modeling are listed below. Some function's headers are displayed.

Cochlear Models

Transmission Line Model (Long Wave Model)

Functions: TlineME, TlineVs, TlinePs
Description:
"Transmission Line" Pressure/Force computation given the distributed impedances for N cochlear sections. All input parametrs are (may be) spectrum vectors. The function TlineVs takes stapes velocity spectrum as input. TlinePs assumes stapes pressure spectrum as input. TLineME assumes a Thevenin's equivalent circuit of external/middle ear in terms of the pressure at eardrum and the ear impedance, both referred to the input of the cochlea. The last section parallel impedance is the helicotrema impedance. The method consists in computing the equivalente line impedances from right-to-left and then use them to compute the pressures/forces at each section.

TLineME:
% function [P,Ps,Zin,Zeq]=TlineME(Pin,Zm,Zs,Zp)
% Computes cochlear pressures and the input impedance of the
% equivalent transmission line cochlear circuit (long wave approx.),
% with middle-ear connected.
% Input: series and parallel distributed impedances:
% series impedance: Zs: Zs(Nf,1) or Zs(Nf,Nx),Nf: no. points in freq.
% parallel impedance: Zp(Nf,Nx), Nx: no. points of x (sections).
% Zm: Middle ear impedance referred to stapes (point/volume)velocity.
% Pin: Input pressure/force referred to stapes.
% Output: Pressures P(1)..P(N); Ps - pressure at stapes
% Input section impedances Zeq; Input coclear impedance, Zin=Zeq(1)
%------------------------------------------------------------------------
% Example with N=5:
%
%      Zm   Ps  Zs1  P1   Zs2 P2  Zs3 P3   Zs4 P4  Zs5   P5
%  +--####--+--@@@@--|--@@@@--|--@@@@--|--@@@@--|--@@@@--+
%  |                 #        #        #        #        #
% (+) Pin        Zp1 #    Zp2 #    Zp3 #    Zp4 #    Zp5 #
%  |                 # |->    #        #        #        #
%  +________.________|_|______|________|________|________+
%          x=0  x=Dx   |    x=2Dx    x=3Dx     x=4Dx    x=5Dx=Lx
%                      |
%                     Zeq(2)
% Ps=P(x=0); P6=P(N)=P(x=L); Zin = Zeq(1)
% Related functions: TLinePs TLineVs TlineME2

Impedance Computation

AllenZ
% function [Zp,Ht]=AllenZ(f,x);
% Computes cochlear partition impedance for Allen's model given
% Nf frequency points in f (column) and Nx sections at positions x.
% Returns also Ht=VC/VBM, the relationship between cilar and
% basilar membrane velocity.
%
% Ref.:J. Allen, "Cochlear Micromechanics-A Physical Model of Transduction"
% JASA 68(6),pp.1660-1669,1980. Also JASA 89(1),pp.287-309, 1991.
% see also Allen.m

Geis93Z
% function [ZG,ZGmod,WQ,M]=Geis93Z(f,x,gamma)
% Computes the partition impedance ZG(f,x) for GEISLER's Model (1993),
% and a related simplified impedance: ZGmod = Za*Ha
% Ref: Hearing Research, 68(1993) 253-262
% Input:
% f - column vector of response's desired frequencies.
% x - row vector of cochlear positions.
% gamma - feedback gain (default: gamma=1.29)
% Output:
% ZG - Partition Impedance of Geisler's model ((Nf x Nx) elements)
% ZGmod - Simplified impedance ZGmod=Za*Ha (Nf x Nx)
% WQ - normalized angular frequencies and quality factors of ZGmod
% WQ=[WZa,WHz,WHp,QZa,QHz,QHp] (6 elements)
% Definitions:
% ZG=ZBM+ZOHC; ZOHC=(Kp/s)/(1+F/Ks)
% ZBM = s*M+R+Kb/s; Kb=M*w0^2; R=d*Kb/w0; ZBM=(M/s)(s^2+s*d*w0+w0^2)
% Kp=Kr*Kc/Ks; Ks=Kr+Kc; Kr=1.17Kb; Kc=0.59Kb;
% Kp=0.392216*Kb = rpb*Kb; Ks=1.76*Kb=rsb*Kb
% Ffb= gamma*Kb*[(wr-s)/(wr+s)]*exp(-s*T); T=n*pi/w0
% With u=jw/w0, ZG/w0 scales, i.e. is independent of x.
% Alternate impedance: ZGmod=Za*Ha
%          M                            (s^2+s*wHz/QHz+wHz^2)
% ZGmod = --- * (s^2+s*wZa/QZa+wZa^2) * ---------------------
%          s                            (s^2+s*wHp/QHp+wHp^2)
%        |<-----------Za------------>| |<---------Ha--------->|
% see also Geis91Z.m


KanisZ
%function [ZKB,WQ,M]=KBoerZ(f,x, gamma)
% Kanis and de Boer Cochlear Model (linear active case)
% Ref.: JASA 94(6), pp. 3199-3206; JASA 96(4), pp. 2156-2165
% Computation of partition impedance ZKB(f,x)
% Input:
% f - column vector of response's desired frequencies
% x - row vector of cochlear positions
% gamma - feedback gain
% Output:
% ZKB - Partition Impedance of Kanis-Boer model ((Nf x Nx) elements)
% WQ - normalized angular frequencies and quality factors of ZKB
% WQ=[WZa,WHz,WHp,QZa,QHz,QHp] (6 elements)
% Note: ZKB scals, i.e. ZKB(s/w0) is independent of x
% ZKB angular freq.s are: WZa*w0, etc.
% Definitions:
% ZBM = (M/s)*(s^2 + s*delta*w0 + w0.^2);
% ZOHC = gamma*c0*(s*w0*(s+w0)/(s^2 + s*dSC*w0 + (sigma*w0)^2);
% ZKB = ZBM-ZOHC ; gamma*M=c0
% w0 = wM*exp(-alfa*x); wM=sqrt(K(x=0)/M)
% ZKB = Za*Ha;
%       M
% Za = --- * (s^2 + s*WZa*w0/QZa + (WZa*w0)^2)
%       s
%
%      (s^2 + s*WHz*w0/QHz + (WHz*w0)^2)
% Ha = ----------------------------------
%      (s^2 + s*WHp*w0/QHp + (WHp*w0)^2)


NKim86Z
% function [Zp, Zpmod, fQ, Mef] = NKim86Z(f, x, gamma)
% Neely and Kim's, 1986 model (JASA 79(5) May 86)
% Computation of partition impedance Zp(f,x)
% Input: % f - column vector for required frequencies
% x - row vector for x positions
% gamma - feedback gain (default: gamma=1)
% Output:
% Zp - Model partition impedance((Nf x Nx) elements)
% Zpmod - Partition impedance as: Zpmod=Za*Ha (Nf x Nx)
% fQ - Parameters of Zpmod: frequencies and quality factors:
% fQ=[fZa;fHz;fHp;QZa;QHz;QHp] (matrix of 6 x Nx elements)
%
%         Mef                         (s^2+s*wHz/QHz+wHz^2)
% Zpmod = ---*(s^2+s*wZa/QZa+wZa^2) * ---------------------
%          s                          (s^2+s*wHp/QHp+wHp^2)
%
%         |<---------Za----------->| |<---------Ha-------->|
%
% Mef = (g/b)*m1 =0.0075 : equivalent mass


Neely93Z
% function [Zp,Hi] = Neely93Z(f, x, gamma)
% Neely's(1993) Coclear Model. Computation of partition impedance Zp(f,x)
% Ref: JASA 94 (1) July 93 pp 137
%
% Input:
% f - column vector for required frequencies
% x - row vector for x positions (must be between 0 and L=2.5cm)
% gamma - feedback gain (default: gamma=1)
% Output:
% Zp - Model Partition SPECIFIC Impedance ((Nf x Nx) elements)
% Hi - IHC gain function

N93ZL
% function Zh=N93Zh(f,gamma);
% Computes helicotrema impedance for no reflections in Neely's model
% (infinite length cochlea)
% Input: f - freq. column vector
% gamma- feedback coefficient (default: 1)
% Output:ZL - helicotrema impedance
%


N93roots
function [wZa,wHz,wHp,aZa,aHz,aHp,wz1,wz2,wp1,wp2]=N93roots(x,gamma)
% Neely's 1993 Model (JASA 94 (1) July 93 pp 137).
% Computation of the ZEROS and POLES of the partition impedance Zp(f,x)
% Input:
% x - row vector for x positions (must be between 0 and L=2.5cm)
% gamma - feedback gain (default: gamma=1)
% Output:
% wZa;wHz;wHp - angular frequencies of quadratic factors
% aZa;aHz;aHp - aZa=wZa/Qza, etc, where Q's are the quality factors
% wz1,wz2,wp1,wp2 - single order zeros and poles of the impedance
%
% Definitions:
%       Mb  (s^2+s*aZa+wZa^2)(s^2+s*aHz+wHz^2)    (s+wz1)(s+wz2)
% Zp = --- * ---------------------------------- * --------------
%       s             (s^2+s*aHp+wHp^2)           (s+wp1)(s+wp2)
%
%NOTE: Zp has 6 zeros and 5 poles. 2 poles and 2 zeros are always real.
%

Example files/Filterbank models

Geis93Ex
% Geisler's 1993 model
% Example
%
% 1. Computs Geisler's 1993 Model Impedance
% 2. Compares impedance with a 4th order rational function impedance
% 3. Computs responses with a Long-Wave model


AllenEx
% Computes cochlear pressure, P, cilia velocity, Vc, with Allen's model.
% see also AllenZ.m
% Ref: S. Puria, J. Allen, "A Parametric Study of Cochlear Input Impedance",
% JASA 89(1), pp. 287-309, Jan. 1991.

KanisEx
%---- Kanis-de Boer's model (linear active case)
% Computes Pressures and Velocities for the Kanis & de Boer' model
% See also: KansiZ.m KanisG.m

NKim86Ex
% Neely and Kim's (1986) cochlear model
% Computes Pressures and Velocities
% Ref: S. Neely, D. Kim, "A Model for Active Elements in Cochlear Biomechanics",
%J. Acoust. Soc. Am. 79(5), pp. 1472-1480, May 1986.

Neely93
% Neely's Model (1993)
% Computation of Responses
% Ref: S. Neely, "A Model of Cochlear Mechanics with Outer Hair Cell Motility",
% J. Acoust. Soc. Am. 94(1), pp. 137-146, July 1993.

NeelyDi
% function [Di,f,x,Zc]=NeelyDi(Nf,fmin,fmax,gamma,Zh);
%
% Computs IHC displacement frequency curves for the Neely's cochlear model
% Returns also the input impedance of the cochlea
% Input: gamma (0.2-1.0) (default 1)
%        Nf - no. of frequency points (default 300)
%        fmin - lowest freq. (default 100 Hz)
%        fmax - highest freq. (default 30000 Hz)
%        Zh - helicotrema impedance (optional)
%
% Output:Di - IHC displacement
%        f - vector of Nf freq. points (log distributed)
%        x - vector of 500 points between 0-L (linear scale)
%        Zc - Cochlear input impedance
%
% Example: [Di,f]=NeelyDi(220); semilogx(f,db(Di(:,1:50:500)))
% see also Neely93Z TLineME Neely93

NeelyCF
% function [f0,x]=NeelyCF(Nx,gamma,Zh)
%
% Computs the characteristic frequencies in the Neely's 1993 model
% as a function of x, the cochlear position
% The characteristic frequencies are defined as the frequencies for which
% the IHC displacement responses peak (Best Frequencies).
%
% Input vars:
% Nx (Number of cochlear sections. Default: 500)
% gamma (default 1) % Zh (optional, the helicotrema impedance (for no reflexions))
%
% NOTE1: Uses 220 frequency values from 100Hz to 30000 Hz.
% NOTE2: Uses polynomial interpolation to smoth the f0 curve
% To see CF´s as a function of x, use: [f0,x]=NeelyCF; plot(x,f0)
% To compare f0 with Liberman function (cat data), use:
% fcat=456*(10 .^(2.1*(1-x/L)) -0.8); plot(x,log10(f0),x,log10(fcat),':')

Kates
% Computs Kates' Auditoy Model freq. responses.
% Ref: J. Kates, "A Time-Domain Digital Cochlear Model",
% IEEE, Signal Proc., Vol.39, No.12, pp. 2573-2592, Dec. 1991
% Note: needs Signal Processing Toolbox (freqz.m)

Rosowski
% function [Zeq,Peq,PC,PS,PT,PEC,PEX,ZiT] = rosowski(f,PPW)
% J. ROSOWSKI's EXTERNAL AND MIDDLE-EAR MODEL
% "Model of External and Middle-Ear Function"
% Ref: "Auditory Computation", SHAR-6, Springer, 1993
% Inputs:
% f - vector with desired frequency points of responses.
% PPW - Plane wave pressure (default=1).
% Outputs:
% PEX - external pressure (at pinna)
% PEC - pressure at ear canal
% PT - pressure at timpani membrane
% PC - pressure at input of cochlea
% Peq - Thevenin pressure;
% Zeq - Thevenin impedance.

CocHum
% function [P,D,x,w0]=CocHum(f)
%
% Cochlear model for the human coclea.
% Uses:
% parametric partition impedance.
% Rosowski external and middle ear model (see rosowski.m)
% Area of cochlear sections as an interpolated polynomial from
% measured data (Weber'49 in Puria and Allen, JASA 89(1),pp.287-309, 1991.)
%
% P: Pressure across cochlear partition
% D: IHC cilia displacement
% x: Coclear distance (500 points from 0 to L=3.5 cm)
% w0: Cochlear map (Liberman function)
%
% Ref:
% "Modelo computacional da coclea humana",
% F. Perdigao, L. Sa', Acustica'98, Lisboa, 1998


CocHumFB
% function H=CocHumFB(W);
%
% Filterbank responses according a cochlear model for the
% human cochlea and for a samplig frequency of 8kHz.
% (see CocHum.m)

SenFB8
% function [HSen, Bpe8, B8, A8]=SenFB8(W);
%
% Computs the freq. responses of fiters in Seneff's model, stage 1
% for a sampling rate of 8kHz. The responses corresponds to a LMS
% fitting of original responses (for 16 kHz) of the first 35 filters
% of the 40 filter model.
% See Slaney's Auditory toolbox, function "SeneffEar.m"
%
% Returns also the filterbank definition for 8kHz as:
% Bpe8: FIR preemphasis filter (order 2)
% B8, A8: 6th order transfer function coefficients.
% NOTE: Chanells 16 to 35 are 5th order.
%
% Example:
%=========
% fs=8000;
% f=logspace(log10(50),log10(fs/2),500)';
% W=2*pi*f/fs;
% H=SenFB8(W);
% semilogx(f,db(H)), pause, semilogx(f,fase(H))

gammaFB
% function [H,BGAMMA,AGAMMA,CF]=gammaFB(W);
%
% Definition of a gammatone filterbank model for
% 8kHz sampling frequency. 35 channels.
% Uses function MakeERBFilters from Slaney's Auditory Toolbox


CocF
% Cochlear (Functional) Model
% Computes an auditory model for 8KHz sampling rate
% Uses 35 sections of cascade filters as well as parallel filters
% Output: Filter coeficients:
% Bm, Am, Bt,At for cascade filters
% Bs As for parallel filters

Inner Hair Cell/Synapse Models

MeddisHC
% function [f,q,w,k]=MeddisHC(s,fs,A,B,h0,f_ini)
% Meddis IHC/Synapse model
%
% s: input signal
% fs: sampling frequency
% A,B,h0: model constants that set fspo, fsat and fiber threshold.
% (see MedConst.m)
%
% If f_ini (mean rate output at t=0) is given,
% use initial conditions that set f(t=0)=f_ini,
% else use initial conditions for sustained null input(f_ini=fspo)


SeneffHC
% function f=SeneffHC(x,fs, A, B, GHW, f_ini)
%
% Implements Seneff's hair cell model.
% Original results with fs=16KHz
% Forward differences are used to approx. derivatives.
%
% See also M. Slaney's Auditory Toolbox.
% There are some differences from function "SeneffEar.m" in that toolbox,
% namely the short-term and rapid adaptation (AGC) filters.
%
% A, B, GHW: parameters that define spontaneous and saturated rate
% f_ini: initial value of f(t) (for filter initial conditions)
% see also: SenConst.m


MartHC
% function f=MartHC(x,fs,A,B,fsat,n)
% Martens-Immerssel IHC/Synapse Model
% x: input signal
% fs: sampling frequency
% fspo - spontaneous firing rate
% fsat - saturation firing rate
% n - power in AGC (default: n=2)


SenConst
% function [A,B,GHW,uf,fmax]=SenConst(Vcenter,fspo,fsat,Vdb);
% Evaluates the Seneff's model constants, A, B and GHW
% in order to obtain a fiber model with spontaneous rate fspo,
% saturated rate fsat and a rate-intensity curve with center
% fcenter=(fsat+fspo)/2 at V=Vcenter.
% The rate-threshold is approx.: Vcenter -25dB. (Vcenter in dB).
% The other constants are not changed in order to have approx.
% the same adaptation time constants.
%
% Vdb is an optional vector with amplitude values for the plot
% of the rate-intensity curve. If "uf" and/or "fmax" are specified,
% then the model is evaluated for Vdb values and returns
% uf, the rate curve, and fmax, the maximum mean rate (on 1ms basis)
% for 1kHz input, 2.5ms rise time, and with the model initially at
% rest conditions.


MedConst
% function [A,B,h0,Vcenter]=MedConst(Th,fspo,fsat,NewConst);
% Computes the Meddis' IHC model constants A, B and h0,
% in order to have a model with absolute threshold Th (dB);
% spontaneous firing rate fspo and saturation firing rate fsat.
%
% The other model constants are left unchanged (default values)
% unless new constants are given in vector NewConst.
% NewConst=[g0,y0,x0,l0,r0]
% default constants: g0=2000; y0=5.05; x0=66.31; l0=2500; r0=6580.
%
% Returns also Vcenter, the sinusoidal amplitude (1kHz) that
% produces a firing rate f=(fsat+fspo)/2.
% Ref: "Implementation Details of a Computational Model of the
% Inner Hair Cell/Auditory-nerve Synapse",
% JASA, 87(4) Apr. 1990, pp.1813-1816.


MarConst
% function [A,B,Vcenter]=MarConst(Th,fspo,fsat,n)
% Computes constants A, B in Martens-Immerseel IHC model,
% in order to have a model with absolute threshold Th (dB);
% spontaneous firing rate fspo and saturation firing rate fsat.
% The other model constants are left unchanged (default values)
% Returns also Vcenter, the sinusoidal amplitude at 1kHz that
% produces a firing rate f=(fsat+fspo)/2.
% n - AGC power factor (default: n=2)
% see also MartHC.m

Utility Functions


PLPfilt
% function [H,E]=PLPfilt(N,Nband,fs,flag_pe)
% Defines Nband filters for use in PLP analysis
% Ref.: H. Hermanski, JASA 87(4), 1990
% N: number of freq. points in the range [0..pi] (pi included! Use N=Nfft/2+1)
% Nband: Num. of effective critical band filters (excluding CFs of 0 and pi)
% fs: sampling frequency
% H: matrix of filter responses (Nband x N)
% E: Preemphasis factors (Nband x 1).
% If "flag_pe" is input and is different from 0,
% apply preemphasis factors to the filters.
% NOTE: for RASTA, preemphasis must be computed later,
% so E, the preemphasis function, is also returned.
% Defaults: Nband=16; fs=10000


MelFilt
% function H = melfilt(Nbank,NFFT,lopass,hipass,fs)
% Creates "Nbanks" mel-frequency triangular filters, H(Nbank,NFFT/2+1)
% First row of H corresponds to the lower CF filter.
% Nbank - number of triangular filters for CFs in ]0..pi[
% NFFT - number of FFT points [0..2*pi[
% lopass - lower cut-off of the first filter
% hipass - upper cut-off of the last filter
% fs - sampling frequency
%
% H may be directly applied to a FFT result. Example:
% Y=fft(X); % with X=X(NFFT,Nframes): 1 column per frame
% Y=abs(Y(1:NFFT/2+1,:)).^2; % power spectrum in [0..pi]
% Z=H*Y; % Power Filter Bank result: Z=Z(Nbank,Nframes)
% % (1 column (parameter vector) per frame)
% NOTE: implementation according to HTK_V2.1, i.e. the filters are
% triangular on a mel scale, not on the linear frequency scale.


LRasta
% function [c,P]=LRasta(x,N,M,Nband,fs,p,Q);
% LOG-RASTA-PLP cepstral analysis of matrix of frames x
% Performs L-RASTA analysis on vector x, ...
% with frames of N samples, every M samples, with Nband filters.
% Ref.: "Rasta-PLP Speech Analysis Technique",
% H. Hermansky et al, ICASSP'92, pp. I-121.
% x: signal vector
% N: frame length
% M: frame rate in samples
% Nband: num. of (effective) critical bands in PLP analysis
% fs: sampling frequency
% c: vector with Q cepstral coefficients (including c0).
% [c,P]=... returns also P, the LRASTA equivalent Power Spectrum
% with 128 points in [0..pi], (liftering effect included).
% NOTE: This function computs cepstrum according to
% Morgan's RASTA_2_2 C-code implementation.
% [www.icsi.berkley.edu/ftp/global/pub/speech/rasta/]
% see also PLP.m, PLPfilt.m, JRasta.m


PLP
% PLP spectral analysis of vector x
% Ref.: H. Hermanski, JASA 87(4), 1990
% P=PLP(x,N,M,Nband,fs) Performs PLP analysis on vector x, ...
% with frames of N samples, every M samples, with Nband critical band filters.
% x: signal vector
% N: frame length
% M: frame rate in samples
% Nband: num. of (effective) critical bands
% fs: sampling frequency
% P: auditory power spectrum (Nband+2)x(Nframes) :
% (endpoints corresponding to 0 and pi are copies from adjacent ones)
%
% Other forms:
% [a,P]=PLP(x,N,M,Nband,fs,p);
% computes vector "a" (p+1)x(nSamples): LP coefficients (order p; a(1)=1).
%
% [c,a,P]=PLP(x,N,M,Nband,fs,p,Q);
% computes vector "c" with Q cepstral coefficients
% (including c0, no liftering).
%--- NOTE: N,M,Nband,p,Q must be integers!


writeHTK
% function writeHTK(fname,X,nSamples,sampPeriod,sampSize,parmKind,SWAP)
% Saves matrix X in a file named fname in HTK format.
% Ref.: The HTK Book, S. Young et al, Entropic Cambridge Research Lab. Ltd.
% sampPeriod - frame sampling period in 100ns units
% paramKind - kind of parameters
% SWAP - if SWAP~=0 : exchange LSBytes with MSBytes (usu. UNIX env.)
% X - matrix with (smapSize/4) rows and nSamples columns
% see also readHTK


ceps2lpc
%function [a,QME]=ceps2lpc(c)
% Conversion from lpc-ceptral coefficients to LPC coefficients.
% Ref: J.Markel & A.Gray, "Linear Prediction of Speech", pp 230.
% Input:
% c : cepstral coefficients c(1,:)..c(Q,:), WITH c(1,:)=log(QME)
% Output:
% a : LPC coefficients in columns, with a(1,:)=1;
% QME : LPC quadratic mean error
%


lpc2ceps
% function c=lpc2ceps(a,Q);
% Conversion from lpc coeficients to lpc-ceptral coeficients.
% Ref: J.Markel & A.Gray, "Linear Prediction of Speech", pp 230
% Inputs
% a : vector or matrix with LPC coefficients in columns, with a(1,:)=1;
% Q : number of desired ceptral coeficients
% Output:
% c : cepstral coefficients c(1,:)..c(Q,:)
% NOTE: c does not include the zero term: c(0)=log(QME)/2
% where QME is the Linear Prediction quadratic mean error.


frames
% function [y,Nf,Nfft]=frames(x,N,M,h);
% Blocks a signal x into Nf frames of length N every M samples.
% Applies a window h (1 column). Default: hamming window.
% Output:
% y: matrix of size (Nfft)x(Nf) (each column is a frame), where:
% Nfft: next power 2 of N: Nfft=2^log2(N);
% Nf: number of frames (or columns): Nf = fix((length(x)-N)/M)+1;


deltas
% function D=deltas(X,N)
% Linear regression of vectors (e.g. cepstral vectors).
% X : matrix with columns as feature vectors per column: X(Nbank,Nframes)
% N : window size: uses values from X(:,n-N) to X(:,n+N) (2N+1) values
% (default: N=2)
% D : delta matrix (with the same size as X)
% D(k)= G*sum(n*X(n+k)), n=-N:N
% G=1/sqrt(sum(n^2)), n=-N:N - (normalizing constant)
% NOTE: the first and last N values are copies from the
% adjacent calculated deltas.


lpc2spec
% function PSD=lpc2spec(a,EMQ,NPSD,fs,range);
% Computes the LPC power spectrum given a matrix, "a", with LP coefficients
% (in columns; a(1,:)=1)
% returns matrix PSD (NPSD x Nf) with LP power spectrum, where:
% NPSD - no. of points in frequency, [0 .. pi] (default: 128)
% fs - sampling freq. (Optional; if specified plots sonogram)
% range - for sonogram plot (default=100dB)
% PSD=|EMQ/A(W)|^2, where:
% A(z)=sum(ak*z^(-k)), k=0..p