From: drowe67 Date: Mon, 16 Mar 2015 23:05:59 +0000 (+0000) Subject: some in progress octave scripts X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=85ee7fd25f3d90fb79fdb13b6d4acafda480f294;p=freetel-svn-tracking.git some in progress octave scripts git-svn-id: https://svn.code.sf.net/p/freetel/code@2077 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/lpcauto.m b/codec2-dev/octave/lpcauto.m new file mode 100644 index 00000000..6895ca23 --- /dev/null +++ b/codec2-dev/octave/lpcauto.m @@ -0,0 +1,116 @@ +function [ar,e,k]=lpcauto(s,p,t) +%LPCAUTO performs autocorrelation LPC analysis [AR,E,K]=(S,P,T) +% Inputs: +% s(ns) is the input signal +% p is the order (default: 12) +% t(nf,3) specifies the frames size details: each row specifies +% up to three values per frame: [len anal skip] where: +% len is the length of the frame (default: length(s)) +% anal is the analysis length (default: len) +% skip is the number of samples to skip at the beginning (default: 0) +% If t contains only one row, it will be used repeatedly +% until there are no more samples left in s. +% +% Outputs: +% ar(nf,p+1) are the AR coefficients with ar(1) = 1 +% e(nf) is the energy in the residual. +% sqrt(e) is often called the 'gain' of the filter. +% k(nf,2) gives the first and last sample of the analysis interval + +% Notes: +% +% (1) The first frame always starts at sample s(1) and the analysis window starts at s(t(1,3)+1). +% (2) The elements of t need not be integers. +% (3) The analysis interval is always multiplied by a hamming window +% (4) As an example, if p=3 and t=[10 5 2], then the illustration below shows +% successive frames labelled a, b, c, ... with capitals for the +% analysis regions. Note that the first frame starts at s(1) +% +% 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 ... +% +% For speech processing, it can be advantageous to restrict the analysis regions +% to time intervals when the glottis is closed. +% +% (5) Frames can overlap: e.g. t=[10 20] will use analysis frames of +% length 20 overlapped by 10 samples. +% (6) For speech processing p should be at least 2*f*l/c where f is the sampling +% frequency, l the vocal tract length and c the speed of sound. For a typical +% male (l=17 cm) this gives f/1000. + +% Copyright (C) Mike Brookes 1997 +% Version: $Id: lpcauto.m 713 2011-10-16 14:45:43Z dmb $ +% +% VOICEBOX is a MATLAB toolbox for speech processing. +% Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 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 2 of the License, or +% (at your option) 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 can obtain a copy of the GNU General Public License from +% http://www.gnu.org/copyleft/gpl.html or by writing to +% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +s = s(:); % make it a column vector +if nargin < 2 p=12; end; +if nargin < 3 t=length(s); end; +%if nargin < 4 w='ham'; end; +[nf,ng]=size(t); +if ng<2 t=[t t]; end; +if ng<3 t=[t zeros(nf,1)]; end; +if nf==1 + nf=floor(1+(length(s)-t(2)-t(3))/t(1)); + tr=0; +else + tr=1; +end; + +ar=zeros(nf,p+1); +ar(:,1)=1; +e=zeros(nf,1); + +t1=1; +it=1; +nw=-1; +zp=zeros(1,p); +r=(0:p); +for jf=1:nf + k(jf,1) = ceil(t1+t(it,3)); + k(jf,2) = ceil(t1+t(it,3)+t(it,2)-1); + cs = (k(jf,1):k(jf,2)).'; + nc = length(cs); + pp=min(p,nc); + dd=s(cs); + if nc~=nw + % possibly we should have a window whose square integral equals unity + ww=hamming(nc); nw=nc; + y=zeros(1,nc+p); + c=(1:nc)'; + end + wd=dd(:).*ww; % windowed data vector + y(1:nc)=wd; % data vector with p appended zeros + z=zeros(nc,pp+1); % data matrix + % was previously z(:)=y(c(:,ones(1,pp+1))+r(ones(nc,1),1:pp+1)); + z(:)=y(repmat(c,1,pp+1)+repmat(r,nc,1)); + rr=wd'*z; + rm=toeplitz(rr(1:pp)); + rk=rank(rm); + if rk + if rk 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 diff --git a/codec2-dev/octave/lspvar.m b/codec2-dev/octave/lspvar.m new file mode 100644 index 00000000..aa4e245a --- /dev/null +++ b/codec2-dev/octave/lspvar.m @@ -0,0 +1,57 @@ +% 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 + + + + + + diff --git a/codec2-dev/octave/save_raw.m b/codec2-dev/octave/save_raw.m new file mode 100644 index 00000000..7f17277e --- /dev/null +++ b/codec2-dev/octave/save_raw.m @@ -0,0 +1,7 @@ +% save_raw.m +% David Rowe 9 Feb 2015 + +function s = save_raw(fn,s) + fs=fopen(fn,"wb"); + fwrite(fs,s,"short"); +endfunction diff --git a/codec2-dev/octave/twomixer.m b/codec2-dev/octave/twomixer.m new file mode 100644 index 00000000..0a718cb5 --- /dev/null +++ b/codec2-dev/octave/twomixer.m @@ -0,0 +1,38 @@ +% twomixer.m +% David Rowe Feb 2015 +% +% test program for combining signals from two mixers + +fc=1024; +Fs=8192; +[b a] = cheby2(6,40,[200 3000]/(Fs/2)); + +w1 = 2*pi*(fc/2)/(Fs/2); +w2 = 2*pi*(3*fc/2)/(Fs/2); +w3 = 2*pi*(5*fc/2)/(Fs/2); +wc = 2*pi*fc/(Fs/2); + +t = 0:(Fs-1); +s1 = cos(w1.*t); +s2 = cos(w2.*t); +s3 = cos(w3.*t); +s = s1 + s2 + s3; + +lo1 = cos(wc*t); +lo2 = cos(2*wc*t); +out1 = filter(b,a,s .* lo1); +out2 = filter(b,a,s .* lo2); + +S = fft(s); + +figure(1) +subplot(211) +plot(abs(S)); +subplot(212) +plot(angle(S)) + +figure(2) +subplot(211) +plot(abs(fft(out1))) +subplot(212) +plot(abs(fft(out2)))