ft_02.cc

这是一个从音频信号里提取特征参量的程序

	    Console::put(output);
	    
	    // output input and output vectors
	    //
	    real_input_vector.debug(L"real_input");
	    freq_vector.debug(L"freq_output");
	  }
	  
	  // compute the inverse complex transform
	  //
	  ft5.setInputLength(freq_vector.length());
	  if (!ft5.compute(inverse_vector, freq_vector)) {
	    Error::handle(name(), L"compute", Error::TEST, __FILE__,
			  __LINE__);
	  }
	  
	  // compute the snr for this operation
	  //
	  // get the sum of the squares for the input
	  //
	  double signal_power = real_input_vector.sumSquare();
	  
	  // get the real part of the inverse vector
	  //
	  inverse_vector.real(real_inverse_vector);
	  
	  // generate the difference signal - for the real transform, we
	  // need to remove the 0.0 imaginary parts - we'll just add them
	  // to the noise component
	  //
	  real_diff_vector.assign(real_input_vector);
	  real_diff_vector.setLength(real_inverse_vector.length(), true);
	  real_diff_vector.sub(real_inverse_vector);
	  double noise_power = (double)real_diff_vector.sumSquare();
	  
	  // compute the snr
	  //
	  Double snr = signal_power / noise_power;
	  snr.log10();
	  snr *= Integral::DB_POW_SCALE_FACTOR;
	  
	  // print the snr
	  //
	  if (level_a >= Integral::DETAILED) {
	    
	    String value;
	    String output;
	    
	    // output the current order
	    //
	    value.assign(snr, L"%g");
	    output.debugStr(name(), L"forward to inverse", L"SNR",
			    value);
	    Console::put(output);
	    
	    // output vectors
	    //
	    inverse_vector.debug(L"inverse_output");
	    
	    Float temp_float;
	    temp_float = signal_power;
	    temp_float.debug(L"signal_power");
	    temp_float = noise_power;
	    temp_float.debug(L"noise_power");
	  }
	  
	  // make sure the snr is sufficient
	  //
	  if (snr < min_float_prec) {
	    real_input_vector.debug(L"real_input");
	    inverse_vector.debug(L"inverse_output");
	    String output, value;
	    value.assign(snr, L"%g");
	    output.debugStr(name(), L"forward to inverse", L"SNR",
			    value);
	    Console::put(output);		
	    Error::handle(name(), L"snr too low", Error::TEST, __FILE__,
			  __LINE__);
	  }

	  // print a separating line
	  //
	  if (level_a >= Integral::DETAILED) {
	    Console::put(L"\n");
	  }
		
	  // output the experiment characteristics
	  //
	  if (level_a >= Integral::DETAILED) {
		
	    String value;
	    String output;
		
	    // output the data precision
	    //
	    value.assign(L"float");
	    output.debugStr(name(), L"forward transform",
			    L"data precision", value);
	    Console::put(output);
		
	    // output the current order
	    //
	    value.assign(curr_desired_length);
	    output.debugStr(name(), L"forward transform",
			    L"transform order", value);
	    Console::put(output);

	    // output the data length
	    //
	    value.assign(curr_length);
	    output.debugStr(name(), L"forward transform",
			    L"data length", value);
	    Console::put(output);

	    // output the forward implementation type
	    //
	    output.debugStr(name(), L"forward transform", L"forward impl",
			    IMPL_MAP.getName(ft5.getImplementation()));
	    Console::put(output);
	  }
	      
	  // compute the forward complex transform
	  //
	  ft5.setAlgorithm(curr_forward_algo);
	  ft5.setImplementation(curr_forward_impl);
	  ft5.setDirection(FourierTransform::FORWARD);
	  ft5.setOrder(curr_desired_length);
	  ft5.setInputLength(curr_desired_length);
	  ft5.setOutputLength(curr_desired_length);	      
	  ft5.setOrder(curr_desired_length);
	  if (!ft5.compute(freq_vector, complex_input_vector)) {
	    Error::handle(name(), L"compute", Error::TEST, __FILE__,
			  __LINE__);
	  }
	      
	  // output the experiment characteristics
	  //
	  if (level_a >= Integral::DETAILED) {
		
	    String value;
	    String output;
		
	    // output the current order
	    //
	    output.debugStr(name(), L"inverse transform", L"inverse impl",
			    IMPL_MAP.getName(ft5.getImplementation()));
	    Console::put(output);

	    // output input and output vectors
	    //
	    complex_input_vector.debug(L"complex_input");
	    freq_vector.debug(L"freq_output");
	  }
	      
	  // compute the inverse complex transform
	  //
	  ft5.setAlgorithm(curr_inverse_algo);
	  ft5.setImplementation(curr_inverse_impl);
	  ft5.setDirection(FourierTransform::INVERSE);
	  ft5.setInputLength(freq_vector.length());
	  if (!ft5.compute(inverse_vector, freq_vector)) {
	    Error::handle(name(), L"compute", Error::TEST, __FILE__,
			  __LINE__);
	  }
	      
	  // compute the snr for this operation
	  //  get the sum of the squares for the input
	  //
	  // signal_power = complex_input_vector.sumSquare();
	  //
	  signal_power = 0;
	  for (long i = 0; i < complex_input_vector.length(); i++) {
	    float temp_float = complex_input_vector(i).mag();
	    signal_power += temp_float * temp_float;
	  }
	      
	  // generate the difference signal - for the complex transform,
	  //
	  complex_diff_vector.assign(complex_input_vector);
	  complex_diff_vector.setLength(inverse_vector.length(),
					true);
	  complex_diff_vector.sub(inverse_vector);

	      // noise_power = complex_diff_vector.sumSquare();
	      //
	  noise_power = 0;
	  for (long i = 0; i < complex_diff_vector.length(); i++) {
	    float temp_float = complex_diff_vector(i).mag();
	    noise_power += temp_float * temp_float;
	  }
	      
	  // compute the snr
	  //
	  snr = signal_power / noise_power;
	  snr.log10();
	  snr *= Integral::DB_POW_SCALE_FACTOR;
	      
	      // print the snr
	      //
	  if (level_a >= Integral::DETAILED) {
		
	    String value;
	    String output;
		
	    // output the current order
	    //
	    value.assign(snr, L"%g");
	    output.debugStr(name(), L"forward to inverse", L"SNR",
			    value);
	    Console::put(output);

	    // output vectors
	    //
	    inverse_vector.debug(L"inverse_output");

	    Float temp_float;
	    temp_float = signal_power;
	    temp_float.debug(L"signal_power");
	    temp_float = noise_power;
	    temp_float.debug(L"noise_power");
	  }

	  // make sure the snr is sufficient
	  //
	  if (snr < min_float_prec) {
	    Error::handle(name(), L"snr too low", Error::TEST, __FILE__,
			  __LINE__);
	  }

	  // print a separating line
	  //
	  if (level_a >= Integral::DETAILED) {
	    Console::put(L"\n");
	  }
	}
	else {
	  if (level_a >= Integral::DETAILED) {
	    String tmp(L"skipping DFT test for this length");
	    Console::put(tmp);
	  }
	}

      } // end loop over inverse transform implementations
    } // end loop over forward transform implementations
  } // end of data length loop
  
  // test FFT algorithm
  //  Case: real spectrum -> complex signal -> complex spectrum
  //
  
  // local variables
  //
  FourierTransform ft6;
  long N = 8;
  VectorFloat real_in(N), real_out, imag_out;
  VectorComplexFloat ft_forward_out, ft_inverse_out;
  real_in.assign(L"1, 2, 6, 8, 3, 7, 5, 4, 4, 5, 7, 3, 8, 6, 2, 1");
  
  // inverse transform
  //
  ft6.set(FourierTransform::FFT, FourierTransform::SPLIT_RADIX,
	 FourierTransform::INVERSE);
    
  ft6.setInputLength(N);

  if (!ft6.compute(ft_inverse_out, real_in)) {
    Error::handle(ft6.name(), L"compute", Error::TEST, __FILE__,
		  __LINE__);
  }
    
  // inverse transform
  //
  ft6.setDirection(FourierTransform::FORWARD);
    
  if (!ft6.compute(ft_forward_out, ft_inverse_out)) {
    Error::handle(ft6.name(), L"compute", Error::TEST, __FILE__,
		  __LINE__);
  }

  // compare inverse output to input
  //
  ft_forward_out.real(real_out);
  ft_forward_out.imag(imag_out);
  real_in.setLength(N);  
  if (!real_out.almostEqual(real_in) || !imag_out.almostEqual(0.0)) {
    imag_out.debug(L"imag output");
    ft_forward_out.debug(L"fft forward result");
    Error::handle(name(), L"fft computation error", Error::TEST, __FILE__,
		  __LINE__);
  }
  

  // test DFT algorithm
  //  Case: real signal -> complex spectrum -> complex signal
  //        input length < output length
  //
  
  // local variables
  //
  FourierTransform ft7;
  N = 7;
  IMPLEMENTATION impl[] = { CONVENTIONAL, TRIGONOMETRIC };
  long num_impl = 2;
  real_in.setCapacity(N + 1);
  real_in.assign(L"1, 2, 6, 8, 3");

  for (int i = 0; i < num_impl; i++) {

    // inverse transform
    //
    ft7.set(FourierTransform::DFT, impl[i], FourierTransform::FORWARD);
      
    ft7.setOutputLength(N);
      
    if (!ft7.compute(ft_forward_out, real_in)) {
      Error::handle(ft7.name(), L"compute", Error::TEST, __FILE__,
		    __LINE__);
    }

    // inverse transform
    //
    ft7.setDirection(FourierTransform::INVERSE);
    ft7.setOutputLength(real_in.length());
    if (!ft7.compute(ft_inverse_out, ft_forward_out)) {
      Error::handle(ft7.name(), L"compute", Error::TEST, __FILE__,
		    __LINE__);
    }

    // compare inverse output to input
    //
    ft_inverse_out.real(real_out);
    ft_inverse_out.imag(imag_out);
      
    if (!real_out.almostEqual(real_in) || !imag_out.almostEqual(0.0)) {
      imag_out.debug(L"imag out");
      ft_inverse_out.debug(L"fft inverse result");
      Error::handle(name(), L"fft computation error", Error::TEST, __FILE__,
		    __LINE__);
    }
  }
  
  // test DCT algorithm
  //

  // local variables
  //
  FourierTransform ft8;    
  N = 16;
  IMPLEMENTATION impl_dct[] = { TYPE_I, TYPE_II, TYPE_III, TYPE_IV };
  num_impl = 4;
  VectorFloat dct_forward_out, dct_inverse_out;
  real_in.setLength(N);  
  for (long i = 0; i < N; i++) {
    real_in(i).cos(2 * Integral::PI * i / N / 4);
  }
    
  for (int i = 0; i < num_impl; i++) {

    // forward transform
    //
    ft8.set(FourierTransform::DCT, impl_dct[i], FourierTransform::FORWARD);
    ft8.setOutputLength(N * 4);

    if (!ft8.compute(dct_forward_out, real_in)) {
      Error::handle(ft8.name(), L"compute", Error::TEST, __FILE__,
		    __LINE__);
    }
      
    // inverse transform
    //
    ft8.setDirection(FourierTransform::INVERSE);
    ft8.setInputLength(N);
      
    if (!ft8.compute(dct_inverse_out, dct_forward_out)) {
      Error::handle(ft8.name(), L"compute", Error::TEST, __FILE__,
		    __LINE__);
    }

    if (level_a >= Integral::DETAILED) {

      ft8.debug(L"DCT");
      real_in.debug(L"real_in");
      dct_forward_out.debug(L"dct_forward_out");
      dct_inverse_out.debug(L"dct_inverse_out");
    }

    // compare inverse output to input
    //
    dct_inverse_out.setLength(N);  
    if (!dct_inverse_out.almostEqual(real_in)) {
      real_in.debug(L"real input");
      dct_inverse_out.debug(L"dct inverse result");
      Error::handle(name(), L"dct computation error", Error::TEST, __FILE__,
		    __LINE__);
    }
  }
  
  // reset indentation
  //
  if (level_a > Integral::NONE) {
    Console::decreaseIndention();
  }

  // --------------------------------------------------------------------
  //
  // 5. print completion message
  //
  // --------------------------------------------------------------------

  // reset indentation
  //
  if (level_a > Integral::NONE) {
    Console::decreaseIndention();
  }
  
  if (level_a > Integral::NONE) {
    String output(L"diagnostics passed for class ");
    output.concat(name());
    output.concat(L"\n");
    Console::put(output);
  }

  // exit gracefully
  //
  return true;
}