// Copyright (C) 2014 egsti // // 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 3 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 should have received a copy of the GNU General Public License // along with this program; if not, see . #include #include #include #define my_i std::complex(0.0, 1.0) DEFUN_DLD (goertzel, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{dft_data} =} goertzel(@var{data}, @var{freq_indices})\n\ computes the discrete Fourier transform of data at continuous frequency\n\ indices freq_indices using the second order Goertzel algorithm as described\n\ by Sysel and Rajmic EURASIP Journal on Advances in Signal Processing 2012:56\n\ \n\ data ........... real or complex input data (vector)\n\ freq_indices ... real frequency indices between 0 and length(data)-1\n\ dft_data ....... in general complex output data\n\ \n\ @end deftypefn") { octave_value_list retval; char fcn[] = "goertzel"; int nargin = args.length(); // Check the input arguments if (nargin < 2) { print_usage (); return retval; } // Read input data dim_vector xdims = args(0).dims(); if (xdims(0)!=1 && xdims(1)!=1) { error("%s: input parameter 1 must be a vector", fcn); } int N = xdims(0)*xdims(1); // signal length ComplexMatrix data = args(0).complex_matrix_value(); // Read frequency indices dim_vector idims = args(1).dims(); if (idims(0)!=1 && idims(1)!=1) { error("%s: index parameter 2 must be a vector", fcn); } int K = idims(0)*idims(1); // number of frequencies to compute Matrix omega = args(1).matrix_value(); // Memory allocation for the output data ComplexMatrix dft_data(idims(0), idims(1)); // Computation via second-order system // loop over the particular frequencies double A=0, B=0; std::complex C = 0, D = 0; std::complex s0 = 0, s1 = 0, s2 = 0; std::complex y0 = 0; octave_idx_type k, i; for (k=0; k