--- /dev/null
+% av_imp.m
+% David Rowe Aug 2012
+% Averages the impulse response samples
+
+function imp = av_imp(imp_filename, period_in_secs, st, en)
+ f = fopen(imp_filename,"rb");
+ s = fread(f, Inf, "short")';
+
+ Fs = 8000;
+ n = period_in_secs * Fs;
+
+ [r c] = size(s);
+
+ imp = zeros(1,n);
+ for i=1:n:c-n
+ imp = imp + s(i:i+n-1);
+ endfor
+
+ % user supplies start and end samples after viweing plot
+
+ if (nargin == 4)
+ imp = imp(st:en);
+ end
+
+ % normalise
+
+ imp /= sqrt(sum(imp .^ 2));
+
+ [h w] = freqz(imp, 1, 4000);
+
+ figure(1);
+ clf;
+ plot(imp);
+
+ figure(2);
+ clf;
+ subplot(211)
+ plot(10*log10(abs(h)))
+ subplot(212)
+ plot(angle(h))
+
+endfunction
+
--- /dev/null
+% cbphase.m
+% David Rowe Aug 2012
+% Used to experiment with critical band phase perception and smoothing
+
+function cbphase
+
+ Wo = 100.0*pi/4000;
+ L = floor(pi/Wo);
+
+ A = zeros(1,L);
+ phi = zeros(1,L);
+
+ % three harmonics in this band
+
+ b = 4; a = b-1; c = b+1;
+
+ % set up phases and mags for 2nd order system (see phasesecord.m)
+
+ wres = b*Wo;
+ phi(a) = 3*pi/4 + wres;
+ phi(b) = pi/2 + wres;
+ phi(c) = pi/4 + wres;
+
+ A(a) = 0.707;
+ A(b) = 1;
+ A(c) = 0.707;
+
+ % add linear component
+
+ phi(1) = pi;
+ phi(2:L) = phi(2:L) + (2:L)*phi(1);
+ phi = phi - 2*pi*(floor(phi/(2*pi)) + 0.5);
+
+ N = 16000;
+ Nplot = 250;
+ s = zeros(1,N);
+
+ for m=a:c
+ s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
+ s = s + s_m;
+ endfor
+
+ figure(2);
+ clf;
+ subplot(211)
+ plot((1:L)*Wo*4000/pi, A,'+');
+ subplot(212)
+ plot((1:L)*Wo*4000/pi, phi,'+');
+
+ %v = A(a)*exp(j*phi(a)) + A(b)*exp(j*phi(b)) + A(c)*exp(j*phi(c));
+ %compass(v,"r")
+ %hold off;
+
+ % est phi1
+
+ diff = phi(b) - phi(a)
+ sumi = sin(diff);
+ sumr = cos(diff);
+ diff = phi(c) - phi(b)
+ sumi += sin(diff);
+ sumr += cos(diff);
+ phi1_ = atan2(sumi, sumr)
+ s_v = cos(Wo*(0:(N-1)) + phi1_);
+
+ figure(1);
+ clf;
+ subplot(211)
+ plot(s(1:Nplot));
+ hold on;
+ plot(s_v(1:Nplot),"r");
+ hold off;
+
+ % build (hopefully) perceptually similar phase
+
+ phi_(a) = a*phi1_;
+ phi_(b) = b*phi1_;
+ phi_(c) = c*phi1_;
+
+ s_ = zeros(1,N);
+
+ for m=a:c
+ s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi_(m));
+ s_ = s_ + s_m;
+ endfor
+
+ subplot(212)
+ plot(s_(1:Nplot));
+
+ gain = 8000;
+ fs=fopen("orig_ph.raw","wb");
+ fwrite(fs,gain*s,"short");
+ fclose(fs);
+ fs=fopen("mod_ph.raw","wb");
+ fwrite(fs,gain*s_,"short");
+ fclose(fs);
+
+endfunction
+
kexc = shift(kexc,sh);
kexc = [kexc(1:sh) zeros(1,512-25) kexc(sh+1:25)];
figure(1)
+clf
plot(kexc)
figure(2)
G = fft(kexc);
+subplot(211)
plot((1:256)*(4000/256),unwrap(angle(G(1:256))))
+subplot(212)
+plot(20*log10(abs(G)))
f=fopen("glottal.c","wt");
fprintf(f,"const float glottal[]={\n");
--- /dev/null
+% lpcpf.m
+% David Rowe Aug 27 2012
+% Experiments with LPC post filtering
+
+function lpcpf(ak_filename, f)
+ aks = load(ak_filename);
+
+ ak = aks(f,:);
+ [tmp p] = size(ak);
+ p -= 1;
+
+ A = freqz(1,ak, 4000);
+ AdB = 20*log10(abs(A));
+
+ gamma = 0.5;
+ gammas = gamma .^ (0:p);
+ W = freqz(ak .* gammas,1, 4000);
+ WdB = 20*log10(abs(W));
+
+ beta = 0.2;
+ R = abs(freqz(ak .* gammas, ak, 4000));
+ %P = (R/max(R)) .^ beta;
+ P = R .^ beta;
+ AP = abs(A) .* P;
+
+ eA = sum(abs(A) .^ 2);
+ eAP = sum(AP .^ 2);
+ gain = sqrt(eA/eAP)
+ AP *= gain;
+
+ PdB = 20*log10(P);
+
+ APdB = 20*log10(AP);
+ 10*log10(sum(AP .^ 2))/10*log10(sum(abs(A) .^ 2))
+
+ figure(1);
+ clf;
+ plot(AdB);
+ hold on;
+ plot(WdB,'g');
+ plot(PdB,'r');
+ plot(APdB,'b.');
+ hold off;
+
+endfunction
+
--- /dev/null
+% phasesecord.m
+% David Rowe Aug 2012
+% Used to experiment with aproximations of phase of 2nd order systems
+
+function phasesecord(w,beta)
+
+ a = [1 -2*cos(w)*beta beta*beta];
+ b = 1;
+
+ [h w1] = freqz(b,a);
+
+ figure(1)
+ subplot(211)
+ plot(abs(h))
+ subplot(212)
+ plot(angle(h))
+
+ % for beta close to 1, we approximate 3 dB points as 1-beta above
+ % and below the resonance freq. Note this fails if w=0 as there is a
+ % double pole. Lets sample the freq response at the 3dB points and
+ % w:
+
+ ws = [w-(1-beta) w w+(1-beta)];
+ [h w1] = freqz(b,a,ws);
+
+ % gain as a fraction of max, should be 3dB. Within 1.3 dB or for w > pi/8,
+ % gets innacurate near w=0 due to 2nd pole
+
+ printf("mag measured...:"); printf("% 4.3f ", abs(h)/max(abs(h)));
+
+ % measured angle, 45 deg from angle at w
+
+ printf("\nangle measured.: "); printf("% 5.3f ", angle(h));
+
+ % Our estimate of angle, (pi+w) is phase at resonance, at lower 3dB
+ % phase is pi/4 ahead, at upper 3B pi/4 behind. -pi/2 is contribution of
+ % other pole at at -w to phase
+
+ ph_lower = (pi+w) + pi/4 - pi/2;
+ ph_res =(pi+w) - pi/2;
+ ph_upper = (pi+w) - pi/4 - pi/2;
+ ph_ests = [ph_lower ph_res ph_upper];
+ ph_ests = ph_ests - 2*pi*(floor(ph_ests/(2*pi)) + 0.5);
+ printf("\nangle estimated:"); printf("% 5.3f ", ph_ests);
+ printf("\n");
+endfunction
+
%
% Plot ampltiude modelling information from dump files.
-function plamp(samname, f)
+function plamp(samname, f, samname2)
+ % switch some stuff off to unclutter display
+
+ plot_lsp = 0;
+ plot_snr = 0;
+ plot_vsnr = 0;
+ plot_sw = 0;
+ plot_pw = 0;
+
sn_name = strcat(samname,"_sn.txt");
Sn = load(sn_name);
snr = load(snr_name);
endif
+ % optional second file, for exploring post filter
+
+ model2q_name = " ";
+ if nargin == 3
+ model2q_name = strcat(samname2,"_qmodel.txt");
+ if file_in_path(".",modelq_name)
+ model2q = load(model2q_name);
+ end
+ end
+
Ew_on = 1;
k = ' ';
do
clf;
% s = [ Sn(2*(f-2)-1,:) Sn(2*(f-2),:) ];
s = [ Sn(2*f-1,:) Sn(2*f,:) ];
- size(s)
+ size(s);
plot(s);
axis([1 length(s) -20000 20000]);
plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;r");
axis([1 4000 -10 80]);
hold on;
-% plot((0:255)*4000/256, Sw(f-2,:),";Sw;");
- plot((0:255)*4000/256, Sw(f,:),";Sw;");
-
- if (file_in_path(".",sw__name))
- plot((0:255)*4000/256, Sw_(f,:),";Sw_;g");
- if (Ew_on == 1)
- plot((0:255)*4000/256, Ew(f,:),";Ew_;r");
- endif
- endif
+ if plot_sw
+ plot((0:255)*4000/256, Sw(f,:),";Sw;");
+ end
if (file_in_path(".",modelq_name))
Amq = modelq(f,3:(L+2));
plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" );
- if (file_in_path(".",pw_name))
+ if (file_in_path(".",pw_name) && plot_pw)
plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;c");
endif
signal = Am * Am';
plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
endif
- if (file_in_path(".",snr_name))
+ if file_in_path(".",model2q_name)
+ Amq2 = model2q(f,3:(L+2));
+ plot((1:L)*Wo*4000/pi, 20*log10(Amq2),";Amq2;m" );
+ end
+
+ if (file_in_path(".",snr_name) && plot_vsnr)
snr_label = sprintf(";Voicing SNR %4.2f dB;",snr(f));
plot(1,1,snr_label);
endif
%plot((1:L)*Wo*4000/pi, 20*log10(orig-synth), phase_err_label);
endif
- if (file_in_path(".",lsp_name))
+ if (file_in_path(".",lsp_name) && plot_lsp)
for l=1:10
plot([lsp(f,l)*4000/pi lsp(f,l)*4000/pi], [60 80], 'r');
endfor
%figure(2);
%endif
- % autocorrelation function to research voicing est
-
- if (file_in_path(".",rk_name))
- figure(3);
- plot(Rk(f,:) / Rk(f,1), ";Rk;");
- hold on;
- p = floor(2*pi/Wo);
- plot([p p ], [0 1], 'r');
- hold off;
- %figure(2);
- endif
-
% interactive menu
printf("\rframe: %d menu: n-next b-back p-png q-quit e-toggle Ew", f);
--- /dev/null
+% Copyright David Rowe 2012
+% This program is distributed under the terms of the GNU General Public License
+% Version 2
+%
+% Plot amplitude modelling information from dump files to test and develop
+% LPC post filter.
+
+function pllpcpf(samname, f)
+
+ % switch some stuff off to unclutter display
+
+ plot_lsp = 0;
+ plot_snr = 0;
+ plot_vsnr = 0;
+ plot_sw = 1;
+ plot_pw = 1;
+ plot_pwb = 1;
+ plot_rw = 1;
+
+ sn_name = strcat(samname,"_sn.txt");
+ Sn = load(sn_name);
+
+ sw_name = strcat(samname,"_sw.txt");
+ Sw = load(sw_name);
+
+ sw__name = strcat(samname,"_sw_.txt");
+ if (file_in_path(".",sw__name))
+ Sw_ = load(sw__name);
+ endif
+
+ model_name = strcat(samname,"_model.txt");
+ model = load(model_name);
+
+ modelq_name = strcat(samname,"_qmodel.txt");
+ if (file_in_path(".",modelq_name))
+ modelq = load(modelq_name);
+ endif
+
+ % Pw (LPC synth filter spectrum) before post filter
+
+ pwb_name = strcat(samname,"_pwb.txt");
+ if (file_in_path(".",pwb_name))
+ Pwb = load(pwb_name);
+ endif
+
+ % Rw (Post filter spectrum)
+
+ rw_name = strcat(samname,"_rw.txt");
+ if (file_in_path(".",rw_name))
+ Rw = load(rw_name);
+ endif
+
+ % Pw (LPC synth filter spectrum) after post filter
+
+ pw_name = strcat(samname,"_pw.txt");
+ if (file_in_path(".",pw_name))
+ Pw = load(pw_name);
+ endif
+
+
+ Ew_on = 1;
+ k = ' ';
+ do
+ figure(1);
+ clf;
+ s = [ Sn(2*f-1,:) Sn(2*f,:) ];
+ size(s);
+ plot(s);
+ axis([1 length(s) -20000 20000]);
+
+ figure(2);
+ Wo = model(f,1);
+ L = model(f,2);
+ Am = model(f,3:(L+2));
+ plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;r");
+ axis([1 4000 -10 80]);
+ hold on;
+ if plot_sw
+ plot((0:255)*4000/256, Sw(f,:),";Sw;");
+ end
+
+ if (file_in_path(".",modelq_name))
+
+ Amq = modelq(f,3:(L+2));
+ plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" );
+
+ if (file_in_path(".",pwb_name) && plot_pwb)
+ plot((0:255)*4000/256, 10*log10(Pwb(f,:)),";Pwb;c");
+ endif
+
+ if (file_in_path(".",rw_name) && plot_rw)
+ plot((0:255)*4000/256, 10*log10(Rw(f,:)),";Rw;c");
+ endif
+
+ if (file_in_path(".",pw_name) && plot_pw)
+ plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;c.");
+ endif
+
+ signal = Am * Am';
+ noise = (Am-Amq) * (Am-Amq)';
+ snr1 = 10*log10(signal/noise);
+ Am_err_label = sprintf(";Am error SNR %4.2f dB;m",snr1);
+ plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
+ endif
+
+
+ hold off;
+
+ % interactive menu
+
+ printf("\rframe: %d menu: n-next b-back p-png q-quit", f);
+ fflush(stdout);
+ k = kbhit();
+ if (k == 'n')
+ f = f + 1;
+ endif
+ if (k == 'b')
+ f = f - 1;
+ endif
+
+ % optional print to PNG
+
+ if (k == 'p')
+ figure(1);
+ pngname = sprintf("%s_%d_sn.png",samname,f);
+ print(pngname, '-dpng', "-S500,500")
+ pngname = sprintf("%s_%d_sn_large.png",samname,f);
+ print(pngname, '-dpng', "-S800,600")
+
+ figure(2);
+ pngname = sprintf("%s_%d_sw.png",samname,f);
+ print(pngname, '-dpng', "-S500,500")
+ pngname = sprintf("%s_%d_sw_large.png",samname,f);
+ print(pngname, '-dpng', "-S1200,800")
+ endif
+
+ until (k == 'q')
+ printf("\n");
+
+endfunction
axis;
if (file_in_path(".", phase_name_))
hold on;
- plot((1:L)*Wo*4000/pi, phase_(f,1:L)*180/pi, "g;phase_;");
+ plot((1:L)*Wo*4000/pi, phase_(f,1:L)*180/pi, "g;phase after;");
grid
hold off;
endif
--- /dev/null
+% Copyright David Rowe 2009
+% This program is distributed under the terms of the GNU General Public License
+% Version 2
+%
+% Plot two sparse phase prediction error text files.
+% Generate data from print_pred_error, print_pred_error_sparse_wo_correction1 etc
+
+function plppe(ppe1_file, ppe2_file, f)
+
+ ppe1 = load(ppe1_file);
+ ppe2 = load(ppe2_file);
+
+ std1 = std(nonzeros(ppe1(:,40:80)));
+ std2 = std(nonzeros(ppe2(:,40:80)));
+
+ printf("std dev for %s is %4.3f\n", ppe1_file, std1);
+ printf("std dev for %s is %4.3f\n", ppe2_file, std2);
+
+ figure(1);
+ clf;
+ subplot(211)
+ hist(nonzeros(ppe1(:,40:80)),20);
+ subplot(212)
+ hist(nonzeros(ppe2(:,40:80)),20);
+
+ k = ' ';
+ do
+ figure(2);
+ clf;
+ subplot(211)
+ L = length(nonzeros(ppe1(f,:)));
+ x = (1:L)*4000/L;
+ std1 = std(nonzeros(ppe1(f,:)));
+ legend = sprintf(";std dev %4.3f;", std1);
+ plot(x, nonzeros(ppe1(f,:)),legend);
+ axis([0 4000 -pi pi]);
+ subplot(212)
+ std2 = std(nonzeros(ppe2(f,:)));
+ legend = sprintf(";std dev %4.3f;", std2);
+ plot(x, nonzeros(ppe2(f,:)),legend);
+ axis([0 4000 -pi pi]);
+
+ % interactive menu
+
+ printf("\rframe: %d menu: n-next b-back p-png q-quit ", f);
+ fflush(stdout);
+ k = kbhit();
+ if (k == 'n')
+ f = f + 1;
+ endif
+ if (k == 'b')
+ f = f - 1;
+ endif
+
+ % optional print to PNG
+
+ if (k == 'p')
+ pngname = sprintf("%s_%d",samname,f);
+ png(pngname);
+ endif
+
+ until (k == 'q')
+ printf("\n");
+
+endfunction
--- /dev/null
+% sd.m
+% David Rowe Aug 2012
+% Plots the spectal distorion between twofiles of LPCs. Used for LSP
+% quantisation tuning.
+
+function sd(raw_filename, dump_file_prefix, f)
+
+ ak1_filename = sprintf("%s_ak.txt", dump_file_prefix);
+ ak2_filename = sprintf("%s_ak_.txt", dump_file_prefix);
+ ak1 = load(ak1_filename);
+ ak2 = load(ak2_filename);
+
+ [ak1_r, ak1_c] = size(ak1);
+ [ak2_r, ak2_c] = size(ak1);
+
+ frames = max([ak1_r ak2_r]);
+ sd = zeros(1,frames);
+ Ndft = 512;
+ A1 = zeros(frames, Ndft);
+ A2 = zeros(frames, Ndft);
+
+ % initial helicopter view of all frames
+
+ for i = 1:frames
+ A1(i,:) = -20*log10(abs(fft(ak1(i,:),Ndft)));
+ A2(i,:) = -20*log10(abs(fft(ak2(i,:),Ndft)));
+ sd(i) = sum((A1(i,:) - A2(i,:)).^2)/Ndft;
+ end
+
+ figure(1);
+ clf;
+ subplot(211)
+ fs=fopen(raw_filename,"rb");
+ s = fread(fs,Inf,"short");
+ plot(s);
+ subplot(212)
+ plot(sd);
+
+ % now enter single step mode so we can analyse each frame
+
+ sn_name = strcat(dump_file_prefix,"_sn.txt");
+ Sn = load(sn_name);
+
+ lsp1_filename = sprintf("%s_lsp.txt", dump_file_prefix);
+ lsp2_filename = sprintf("%s_lsp_.txt", dump_file_prefix);
+ lsp1 = load(lsp1_filename);
+ lsp2 = load(lsp2_filename);
+
+ k = ' ';
+ do
+ figure(2);
+ clf;
+ subplot(211)
+ s = [ Sn(2*f-1,:) Sn(2*f,:) ];
+ size(s);
+ plot(s);
+ axis([1 length(s) -20000 20000]);
+
+ subplot(212);
+ plot((1:Ndft/2)*4000/(Ndft/2), A1(f,1:(Ndft/2)),";A1;r");
+ axis([1 4000 -20 40]);
+ hold on;
+ plot((1:Ndft/2)*4000/(Ndft/2), A2(f,1:(Ndft/2)),";A2;");
+
+ for l=1:10
+ plot([lsp1(f,l)*4000/pi lsp1(f,l)*4000/pi], [0 -10], 'r');
+ plot([lsp2(f,l)*4000/pi lsp2(f,l)*4000/pi], [-10 -20], 'b');
+ endfor
+ plot(0,0,';lsp1;r');
+ plot(0,0,';lsp2;b');
+ sd_str = sprintf(";sd %3.2f dB*dB;", sd(f));
+ plot(0,0,sd_str);
+
+ hold off;
+
+ % interactive menu
+
+ printf("\rframe: %d menu: n-next b-back q-quit", f);
+ fflush(stdout);
+ k = kbhit();
+ if (k == 'n')
+ f = f + 1;
+ endif
+ if (k == 'b')
+ f = f - 1;
+ endif
+
+ until (k == 'q')
+ printf("\n");
+
+endfunction
+
--- /dev/null
+% twotone.m
+% David Rowe Aug 2012
+% Used to experiment with combining phase of two tones
+
+function cbphase
+
+ Wo = 100.0*pi/4000;
+ L = floor(pi/Wo);
+ phi = zeros(1,L);
+
+ % two harmonics
+
+ a = 20; b = 21;
+
+ % set up phases to whatever
+
+ phi(a) = -pi;
+ phi(b) = -pi/2;
+
+ % synthesis the two-tone signal
+
+ N = 16000;
+ Nplot = 250;
+ s = zeros(1,N);
+
+ for m=a:b
+ s_m = cos(m*Wo*(0:(N-1)) + phi(m));
+ s = s + s_m;
+ endfor
+
+ % now our theory says that this signal should be the same perceptually
+
+ phi_(a) = (phi(a) - phi(b))/2;
+ phi_(b) = (phi(b) - phi(a))/2;
+
+ s_ = zeros(1,N);
+ for m=a:b
+ s_m = cos(m*Wo*(0:(N-1)) + phi_(m));
+ s_ = s_ + s_m;
+ endfor
+
+ % plot them and see if envelope has the same phase, but "carriers"
+ % have different phase
+
+ figure(1);
+ clf;
+ subplot(211);
+ plot(s(1:Nplot));
+ subplot(212);
+ plot(s_(1:Nplot),'r');
+endfunction
+
--- /dev/null
+% twotone1.m
+% David Rowe Aug 17 2012
+%
+% Used to experiment with combining phase of two tones. This version
+% sets up a complete synthetic speech signal then tries to combine the
+% phase of high frequency tones. Lets see if we can do this and keep perceptual
+% nature of signal the same.
+
+function twotone1
+
+ % hts1a frame 47
+
+ Wo = 0.093168;
+ L = 33;
+ %A = [69.626907 460.218536 839.677429 2577.498047 972.647888 712.755066 489.048553 364.830536 409.230652 371.767487 489.112854 893.127014 2447.596680 752.878113 475.720520 234.452271 248.161606 232.171051 202.669891 323.914490 678.749451 362.958038 211.652512 170.764435 148.631790 169.261673 272.254150 176.872375 67.344391 99.022301 60.812035 34.319073 14.864757];
+ A = zeros(1,L)*100;
+ phi = [1.560274 1.508063 -1.565184 1.289117 -2.547365 1.412528 -1.303992 3.121130 1.087573 -1.158161 -2.928007 0.995093 -2.614023 0.246136 -2.267406 2.143802 -0.273431 -2.266897 1.685171 -0.668712 2.699722 -1.151891 2.406379 -0.046192 -2.718611 0.761067 -2.305014 0.133172 -1.428978 1.492630 -1.668385 1.539734 -1.336615];
+ %phi = zeros(1,L);
+ st = floor(L/2);
+ %st = 1;
+
+ A(st:st+5) = 1000;
+
+ % now set up phase of signal with phase of upper frequency harmonic
+ % pairs combined
+
+ phi_ = phi;
+ for m=floor(L/2):2:L
+ phi_(m) = (phi(m) - phi(m+1))/2;
+ phi_(m+1) = (phi(m+1) - phi(m))/2;
+ %phi_(m+1) = 0;
+ end
+
+ % synthesise the signals
+
+ N = 16000;
+ Nplot = 250;
+
+ s = zeros(1,N);
+ for m=st:L
+ s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
+ s = s + s_m;
+ endfor
+
+ s_ = zeros(1,N);
+ for m=st:L
+ s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi_(m));
+ s_ = s_ + s_m;
+ endfor
+
+ % plot them, expect to see similar time domain waveforms
+
+ figure(1);
+ clf;
+ subplot(211);
+ plot(s(1:Nplot));
+ subplot(212);
+ plot(s_(1:Nplot),'r');
+
+ figure(2);
+ clf;
+ subplot(211);
+ plot(s(1:Nplot)-s_(1:Nplot));
+
+ % save to disk
+
+ gain = 1;
+ fs=fopen("twotone1_orig.raw","wb");
+ fwrite(fs,gain*s,"short");
+ fclose(fs);
+ fs=fopen("twotone1_comb.raw","wb");
+ fwrite(fs,gain*s_,"short");
+ fclose(fs);
+
+endfunction
+
projects / freetel-svn-tracking.git / commitdiff