--- /dev/null
+function [ar,e,k]=lpcauto(s,p,t)\r
+%LPCAUTO performs autocorrelation LPC analysis [AR,E,K]=(S,P,T)\r
+% Inputs:\r
+% s(ns) is the input signal\r
+% p is the order (default: 12)\r
+% t(nf,3) specifies the frames size details: each row specifies\r
+% up to three values per frame: [len anal skip] where:\r
+% len is the length of the frame (default: length(s))\r
+% anal is the analysis length (default: len)\r
+% skip is the number of samples to skip at the beginning (default: 0)\r
+% If t contains only one row, it will be used repeatedly\r
+% until there are no more samples left in s.\r
+%\r
+% Outputs:\r
+% ar(nf,p+1) are the AR coefficients with ar(1) = 1\r
+% e(nf) is the energy in the residual.\r
+% sqrt(e) is often called the 'gain' of the filter.\r
+% k(nf,2) gives the first and last sample of the analysis interval\r
+\r
+% Notes:\r
+%\r
+% (1) The first frame always starts at sample s(1) and the analysis window starts at s(t(1,3)+1).\r
+% (2) The elements of t need not be integers.\r
+% (3) The analysis interval is always multiplied by a hamming window\r
+% (4) As an example, if p=3 and t=[10 5 2], then the illustration below shows\r
+% successive frames labelled a, b, c, ... with capitals for the\r
+% analysis regions. Note that the first frame starts at s(1)\r
+%\r
+% a a A A A A A a a a b b B B B B B b b b c c C C C C C c c c d ...\r
+%\r
+% For speech processing, it can be advantageous to restrict the analysis regions\r
+% to time intervals when the glottis is closed.\r
+%\r
+% (5) Frames can overlap: e.g. t=[10 20] will use analysis frames of\r
+% length 20 overlapped by 10 samples.\r
+% (6) For speech processing p should be at least 2*f*l/c where f is the sampling\r
+% frequency, l the vocal tract length and c the speed of sound. For a typical\r
+% male (l=17 cm) this gives f/1000.\r
+\r
+% Copyright (C) Mike Brookes 1997\r
+% Version: $Id: lpcauto.m 713 2011-10-16 14:45:43Z dmb $\r
+%\r
+% VOICEBOX is a MATLAB toolbox for speech processing.\r
+% Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html\r
+%\r
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
+% This program is free software; you can redistribute it and/or modify\r
+% it under the terms of the GNU General Public License as published by\r
+% the Free Software Foundation; either version 2 of the License, or\r
+% (at your option) any later version.\r
+%\r
+% This program is distributed in the hope that it will be useful,\r
+% but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\r
+% GNU General Public License for more details.\r
+%\r
+% You can obtain a copy of the GNU General Public License from\r
+% http://www.gnu.org/copyleft/gpl.html or by writing to\r
+% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.\r
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
+\r
+s = s(:); % make it a column vector\r
+if nargin < 2 p=12; end;\r
+if nargin < 3 t=length(s); end;\r
+%if nargin < 4 w='ham'; end;\r
+[nf,ng]=size(t);\r
+if ng<2 t=[t t]; end;\r
+if ng<3 t=[t zeros(nf,1)]; end;\r
+if nf==1\r
+ nf=floor(1+(length(s)-t(2)-t(3))/t(1));\r
+ tr=0;\r
+else\r
+ tr=1;\r
+end;\r
+\r
+ar=zeros(nf,p+1);\r
+ar(:,1)=1;\r
+e=zeros(nf,1);\r
+\r
+t1=1;\r
+it=1;\r
+nw=-1;\r
+zp=zeros(1,p);\r
+r=(0:p);\r
+for jf=1:nf\r
+ k(jf,1) = ceil(t1+t(it,3));\r
+ k(jf,2) = ceil(t1+t(it,3)+t(it,2)-1);\r
+ cs = (k(jf,1):k(jf,2)).';\r
+ nc = length(cs);\r
+ pp=min(p,nc);\r
+ dd=s(cs);\r
+ if nc~=nw\r
+ % possibly we should have a window whose square integral equals unity\r
+ ww=hamming(nc); nw=nc;\r
+ y=zeros(1,nc+p);\r
+ c=(1:nc)';\r
+ end\r
+ wd=dd(:).*ww; % windowed data vector\r
+ y(1:nc)=wd; % data vector with p appended zeros\r
+ z=zeros(nc,pp+1); % data matrix\r
+ % was previously z(:)=y(c(:,ones(1,pp+1))+r(ones(nc,1),1:pp+1));\r
+ z(:)=y(repmat(c,1,pp+1)+repmat(r,nc,1));\r
+ rr=wd'*z;\r
+ rm=toeplitz(rr(1:pp));\r
+ rk=rank(rm);\r
+ if rk\r
+ if rk<pp\r
+ rm=rm(1:rk,1:rk);\r
+ end\r
+ ar(jf,2:rk+1)=-rr(2:rk+1)/rm;\r
+ end\r
+ e(jf)=rr*ar(jf,1:pp+1)';\r
+ t1=t1+t(it,1);\r
+ it=it+tr;\r
+end\r
+\r
--- /dev/null
+%
+% Given input AR filter coefficients, this function calculates
+% the corresponding Line spectral pair (LSP) parameters.
+%
+% Author : Madhukar Budagavi
+%
+
+1;
+
+function w=atolsp(a)
+
+ [m,n] = size(a);
+ if m > n
+ a = a';
+ end
+ m = length(a)-1;
+ pz = [a,0] - [0 fliplr(a)];
+ qz = [a,0] + [0 fliplr(a)];
+ angpz = sort(angle(roots(pz)));
+ angpz = angpz(find(angpz > 0));
+ angqz = sort(angle(roots(qz)));
+ angqz = angqz(find(angqz > 0 & angqz < pi));
+ w = zeros(1,m);
+ j = 1:(m/2);
+ w(2*j) = angpz';
+ w(2*j-1) = angqz';
+endfunction
+
+%
+% Given input Line spectral pair (LSP) parameters of an AR filter,
+% this function calculates the corresponding AR filter parameters
+%
+% Author : Madhukar Budagavi
+%
+
+function a=lsptoa(w)
+
+ [m,n] = size(w);
+ if m > n
+ w = w';
+ end
+ m = length(w);
+ pz = 1;
+ qz = 1;
+ for j = 1:(m/2),
+ pz = conv(pz,[1 -2*cos(w(2*j)) 1]);
+ qz = conv(qz,[1 -2*cos(w(2*j-1)) 1]);
+ end
+ pz = conv([1 -1],pz);
+ qz = conv([1 1],qz);
+ a = (pz+qz)/2;
+ a = a(1:(m+1));
+
+endfunction
--- /dev/null
+% lspvar.m
+% David Rowe 8 Feb 2015
+%
+% Experiments in reducing LSP variance in order to make quantisation easier
+
+lsp;
+close all;
+graphics_toolkit ("gnuplot");
+
+function v = det_var(s)
+ p = 10;
+ ak=lpcauto(s,p,[160 80]);
+ [rows cols] = size(ak);
+ w = zeros(rows,p);
+ for r=1:rows
+ w(r,:) = atolsp(ak(r,:));
+ end
+ v = sum(var(w));
+ figure;
+ plot(w)
+end
+
+s=load_raw("../raw/cq_ref.raw");
+det_var(s)
+
+if 1
+ Fs = 8000;
+ [b a] = cheby2(6,40,[300 2600]/(Fs/2));
+ s1 = filter(b,a,s);
+ det_var(s1)
+end
+
+% + Equalise to remove spectral slope, this reduces variance of spectral
+% envelope
+% + Can probably get rid of energy for many vectors too, except flat ones?
+% + Or have some sort of correlation of gain with peakiness of vectors?
+% + Once flattened try slicing out low energy harmonics, like post filter. How
+% does this sound? Any artefacts of harmonics coming and going near
+% thresholds?
+% + BPF as well, it's all abt minimal information.
+% + What's left? Look at it. How can it be represented? Look at tol distort,
+% what quantiser will deliver that? Just position of clipped peaks?
+% + fall back HF mode, no energy, voicing, just spectral information. Diveristy
+% compared to FEC? Test with simulation of ideal code. Aux carrier with
+% extra information to get 1300-ish quality. Coherent, pilot assisted demod.
+% Low PAPR as extra carriers at lower power level. Essentially for free in
+% terms of TX power.
+% + What is a viable plan to get this going in a mon month or so?
+% + Octave equlaisation of input files
+% + take model files, simulate quantisation, back to C for synth?
+% + Coherent FDM modem port to C
+
+
+
+
+
+
projects / freetel-svn-tracking.git / commitdiff