function helix_values = complex_helix(duration, frequency, samplerate);

  timescale = (0:ceil(duration*samplerate)-1) / samplerate;

  helix_values = exp(%i*2*%pi*frequency*timescale);

//end complex_helix

function helix_values = complex_glisshelix(duration, frequency, samplerate, glissrate);

  lasttick = ceil(duration*samplerate)-1;

  timescale = (0:lasttick);

  timescale = (timescale + glissrate * (timescale-lasttick/2)^2) / samplerate;

  helix_values = exp(%i*2*%pi*frequency*timescale);

//end complex_glisshelix

function argout = arg(x);

  argout = atan(imag(x), real(x));

//end arg

function plot_signal_ri(sigin, samprate);

  box_scale = 1.2;
  sigsize = size(sigin, 2);

  axis = (0:sigsize-1) / samprate;
  sigamp = max(max(abs(real(sigin))), max(abs(imag(sigin))));
  plot_box = [axis(1) -sigamp*box_scale axis(sigsize-1) sigamp*box_scale];

  plot2d(axis, real(sigin), 2, "010", " ", plot_box);

  plot2d(axis, imag(sigin), 5, "011", " ", plot_box);

//end plot_signal_ri

function plot_spectrum(specin, min_freq, max_freq);

  box_scale = 1.2;
  plot_gap = 0.2;
  specsize = size(specin,2);

  axis = (0:specsize-1)*(max_freq-min_freq)/(specsize-1)+min_freq;
  specamp = max(abs(specin));
  plot_box = [min_freq -(1+plot_gap)*specamp*box_scale ...
              max_freq specamp*box_scale];

  plot2d(axis, zeros(specin), 5, "011", " ", plot_box);
  plot2d(axis, -plot_gap*specamp*ones(specin), 5, "010", " ", plot_box);
  plot2d(axis, -(0.5+plot_gap)*specamp*ones(specin), 5, "010", " ", plot_box);
  plot2d(axis, -(1+plot_gap)*specamp*ones(specin), 5, "010", " ", plot_box);

  plot2d(axis, abs(specin), 2, "010", " ", plot_box);
  plot2d(axis, (arg(specin)-%pi)*specamp/(2*%pi)-plot_gap*specamp, 15, ...
         "010", " ", plot_box);

//end plot_spectrum

function plot_spec_path(specin);

  box_scale = 1.2;

  max_coord = max(max(abs(real(specin))), max(abs(imag(specin))));

  plot_box = [-max_coord -max_coord max_coord max_coord] * box_scale;

  plot2d([0 0], [-max_coord max_coord], 5, "011", " ", plot_box);
  plot2d([-max_coord max_coord], [0 0], 5, "010", " ", plot_box);

  plot2d(real(specin), imag(specin), 2, "010", " ", plot_box);
  plot2d(real(specin), imag(specin), -15, "010", " ", plot_box);

//end plot_spec_path

function specout = spectrum(sigin, number_bins);

  zeropad = max(number_bins-size(sigin,2), 0);

  specout = fft([sigin zeros(1, zeropad)], -1);

//end spectrum

function plot_zc_trim_spectrum(sigin, freqres, ...
                               samprate, minfreq, maxfreq);

  sigsize = size(sigin, 2);

  number_freq_bins = ceil(samprate / freqres);

  true_freqres = samprate / number_freq_bins;

  zero_freq_bin = ceil(number_freq_bins/2);

  min_out_bin = floor(minfreq / true_freqres) + zero_freq_bin;

  max_out_bin = ceil(maxfreq / true_freqres) + zero_freq_bin;

  true_minfreq = (min_out_bin - zero_freq_bin) * true_freqres;

  true_maxfreq = (max_out_bin - zero_freq_bin) * true_freqres;

  specplot = mtlb_fftshift(spectrum(sigin, number_freq_bins));

  specplot = specplot(min_out_bin:max_out_bin);

  plot_spectrum(specplot, true_minfreq, true_maxfreq);

//end plot_zc_trim_spectrum

function specout = shift_phase(specin, shift_samp);

  spec_length = size(specin,2);

  t = (0:spec_length-1);

  specout = specin .* exp(-%i*2*%pi*(shift_samp/spec_length)*t);

//end shift_phase

function ft_helix_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = complex_helix(duration, frequency, samplerate);

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_helix_440_demo

function ft_helix_0_demo;

  samplerate = 44100;
  frequency = 0;
  duration = 1/110;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = complex_helix(duration, frequency, samplerate);

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_helix_0_demo

function ft_cos_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = real(complex_helix(duration, frequency, samplerate));

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_cos_440_demo

function ft_sin_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = imag(complex_helix(duration, frequency, samplerate));

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_sin_440_demo

function ft_isin_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = %i * imag(complex_helix(duration, frequency, samplerate));

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_isin_440_demo

function strout = stri(compno);

  strout = real(compno) + %i * 2 * imag(compno);

//end stri

function ft_ehelix_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = stri(complex_helix(duration, frequency, samplerate));

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_ehelix_440_demo

function ft_helix_440_990_demo;

  samplerate = 44100;
  fundfreq = 440;
  partial2 = 990;
  duration = 4/fundfreq;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = complex_helix(duration, fundfreq, samplerate) ...
           + complex_helix(duration, partial2, samplerate);

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_helix_440_990_demo

function ft_swellhelix_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = complex_helix(duration, frequency, samplerate) ...
                        .* (0:ceil(duration*samplerate)-1) / samplerate;

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_swellhelix_440_demo

function ft_glisshelix_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  grate = 0.001;
  signal = complex_glisshelix(duration, frequency, samplerate, grate);

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_glisshelix_440_demo

function ft_helix_pulse_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 2;
  minfreq = 300;
  maxfreq = 500;
  freqres = 1 / duration;
  signal = loadwave('helix_spulse_440.wav');

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 2);
  xbasc();
  xselect();

  plot_signal_ri(signal(42000:46000), samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_helix_pulse_440_demo

function ft_gnoise_440_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 2;
  minfreq = 300;
  maxfreq = 500;
  freqres = 1 / duration;
  signal = loadwave('gnoise_440.wav');

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 2);
  xbasc();
  xselect();

  plot_signal_ri(signal(42000:46000), samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_gnoise_440_demo

function ft_helix_440_if_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 1.5/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / duration;
  signal = complex_helix(duration, frequency, samplerate);

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

//end ft_helix_440_if_demo

function ft_helix_440_hfr_demo;

  samplerate = 44100;
  frequency = 440;
  duration = 4/frequency;
  minfreq = -5000;
  maxfreq = 5000;
  freqres = 1 / (10 * duration);
  signal = complex_helix(duration, frequency, samplerate);

  xset("window", 0);
  xbasc();
  xselect();

  plot_signal_ri(signal, samplerate);

  xset("window", 1);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, minfreq, maxfreq);

  xset("window", 2);
  xbasc();
  xselect();

  plot_zc_trim_spectrum(signal, freqres, samplerate, frequency-500, frequency+500);

//end ft_helix_440_hfr_demo

function sincout = mysinc(sincin);

  insize = size(sincin, 2);

  sincout = zeros(sincin);

  for i = 1:insize,

    if sincin(i) == 0 then
      sincout(i) = 1;
    else
      sincout(i) = sin(%pi*sincin(i))/(%pi*sincin(i));
    end,

  end

//end mysinc

function sinc_demo;

  xset("window", 0);
  xbasc();
  xselect();

  minvalue = -10;
  maxvalue = 10;
  stepvalue = 0.05;
  domain = (minvalue/stepvalue:maxvalue/stepvalue)*stepvalue;

  box_scale = 1.2;
  plot_box = [minvalue -1 maxvalue box_scale];

  plot2d(domain, zeros(domain), 5, "011", " ", plot_box);
  plot2d(domain, mysinc(domain), 2, "010", " ", plot_box);

//    plot(sincd(1, 1));

//end sinc_demo
