safiire/lowpass.cpp

## lowpass.cpp
#include <iostream>
#include <complex>

#define PI 3.141592653589793
#define SAMPLE_RATE 64

using namespace std;

//  This lowpass is a weighted moving average
//  If you change the + to - it will become a highpass filter
void lowpass(float a, float b, float *input, float *output, int size){

  output[0] = a * input[0];  // avoid the negative index for n = 0

  for(int n = 1; n < size; n++){
    output[n] = a * input[n] + b * input[n - 1];
  }
}


//  Print out the values of an array.
void print_buffer(float *buffer, int size){
  cout << "[";
  for(int n = 0; n < size; n++){
    cout << buffer[n] << ", ";
  }
  cout << "]\n\n";
}


//  Crappy O(n^2) discrete fourier transform code
//  ∑ e^jΩ * input[n]
void dft(float *input, std::complex<float> *output, int size){
  float root = 2 * 3.141592653589793 * 1 / (float)size;
  for(int bin = 0; bin < size; bin++){
    complex<float> c = polar(1.f, root * bin);
    complex<float> sum = complex<float>(0, 0);
    for(int n = 0; n < size; n++){
      sum += input[n] * std::pow(c, n);
    }
    output[bin] = sum;
  }
}


//  This is some kind of attempt to display a frequency response
void print_frequency_response_graph(float *input, int size){
  const char *line = " -------------------------------NN-------------------------------";

  int quantized[size];
  for(int n = 0; n < size; n++){
    quantized[n] = (int)(input[n] * 10);
  }

  cout << line << endl;
  for(int magnitude = 10; magnitude >= 0; magnitude--){
    cout << "|";
    for(int n = 0; n < size; n++){
      if(quantized[n] == magnitude){
        cout << "*";
      }else{
        cout << " ";
      }
    }
    cout << "|" << endl;
  }
  cout << line << endl;
}

//  This tries to print an audio buffer
void print_audio_buffer(float *input, int size){
  const char *line = " ----------------------------------------------------------------";

  int quantized[size];
  for(int n = 0; n < size; n++){
    quantized[n] = (int)(input[n] * 10);
  }

  cout << line << endl;
  for(int magnitude = 10; magnitude >= -10; magnitude--){
    cout << (magnitude == 0 ? "0" : "|");
    for(int n = 0; n < size; n++){
      if(quantized[n] == magnitude){
        cout << "*";
      }else{
        cout << " ";
      }
    }
    cout << (magnitude == 0 ? "0" : "|") << endl;
  }
  cout << line << endl;
}


int main(void){

  float size = SAMPLE_RATE;

  //  Change these:
  float frequency_hz = 25.f;    // This is the frequency of the input cosine wave

  //  Here are the coefficients
  //  Make a + b = 1 so the filter has no overall gain
  float b = 0.4f;
  float a = 1.f - b;


  //  Make a buffer holding 1 second of audio at the sample rate
  //  and fill it up a 25hz sine wave
  float input[SAMPLE_RATE];
  float w = 2 * PI * (frequency_hz / (float)SAMPLE_RATE); // make this radian frequency

  for(int n = 0; n < size; n++){
    input[n] = cos(w * n);
  }

  //  Buffer to hold the output of the lowpass filter
  float low[SAMPLE_RATE] = {0};

  //  Buffer to hold an impulse, ie [1, 0, 0, 0, 0, 0, ...]
  float impulse[SAMPLE_RATE] = {0};
  impulse[0] = 1.f;

  //  And a buffer to hold the filter's reponse to an impulse
  float impulse_response[SAMPLE_RATE] = {0};


  //  Display the input buffer
  cout << "Input: " << endl;
  print_buffer(input, size);
  print_audio_buffer(input, size);

  //  Low pass filter the input
  lowpass(a, b, input, low, size);

  //  Display the lowpassed output
  cout << "Lowpass Output: " << endl;
  print_buffer(low, size);
  print_audio_buffer(low, size);

  //  Figure out the impulse response of this filter
  //  All we do is send an impulse into the filter and see what comes out
  lowpass(a, b, impulse, impulse_response, size);

  //  Display the impulse response
  cout << "Impulse Response" << endl;
  print_buffer(impulse_response, size);
  print_audio_buffer(impulse_response, size);

  //  A complex valued buffer to hold the filter's frequency response
  std::complex<float> frequency_response[SAMPLE_RATE];

  //  Make a buffer to hold the magnitude reponse
  float magnitude_response[SAMPLE_RATE];

  //  Magically, the fourier transform of the impulse reponse IS
  //  the frequency response
  dft(impulse_response, frequency_response, size);

  //  The frequency response contains the magnitude response if you
  //  take the absolute value of each complex value
  //  This is what we see as a graph on graphic EQs :)
  for(int n = 0; n < size; n++){
    magnitude_response[n] = abs(frequency_response[n]);
  }

  //  Display the magnitude reponse.  But remember that frequency
  //  spectrums of real signals are mirror images about the nyquist (center)
  //  frequency.  Generally we only see the first half of this in programs
  //  so a low pass filter is going to look like a dip down and then back up
  //  in the mirrored half.
  cout << "Magnitude Response" << endl;
  print_buffer(magnitude_response, size);
  print_frequency_response_graph(magnitude_response, size);

  return 0;
}

/* Output

Input:
[1, -0.77301, 0.19509, 0.471397, -0.923879, 0.95694, -0.55557, -0.0980161, 0.707107, -0.995185, 0.83147, -0.290285, -0.382684, 0.881921, -0.980786, 0.634394, -7.35157e-07, -0.634393, 0.980785, -0.881922, 0.382685, 0.290284, -0.831469, 0.995185, -0.707106, 0.0980195, 0.555569, -0.956939, 0.923881, -0.4714, -0.195088, 0.773009, -1, 0.773011, -0.19509, -0.471397, 0.92388, -0.95694, 0.555574, 0.0980127, -0.707104, 0.995184, -0.831471, 0.290286, 0.382682, -0.881921, 0.980785, -0.634393, -1.60923e-06, 0.634389, -0.980784, 0.881923, -0.382687, -0.290282, 0.831464, -0.995185, 0.707113, -0.0980171, -0.555564, 0.956941, -0.923882, 0.471394, 0.195086, -0.773013, ]
 ----------------------------------------------------------------
|*                                                               |
|     *            *    *    *       *    *    *            *    |
|          *  *                                     *  *         |
|        *                      * *                      *       |
|               *                                 *              |
|                          *           *                         |
|   *                                                         *  |
|                    *                       *                   |
|                     *                     *                    |
|  *                                                           * |
0       *        *        *             *        *        *      0
|                              *   *                             |
|           *                                         *          |
|            *                                       *           |
|                             *     *                            |
|      *                                                   *     |
|                 *                             *                |
| *                      *               *                      *|
|                   *  *                   *  *                  |
|    *    *    *            *         *            *    *    *   |
|                                *                               |
 ----------------------------------------------------------------

Lowpass Output:
[0.6, -0.0638063, -0.19215, 0.360874, -0.365769, 0.204613, 0.0494342, -0.281038, 0.385057, -0.314268, 0.100808, 0.158417, -0.345724, 0.376079, -0.235703, -0.0116777, 0.253757, -0.380636, 0.334714, -0.136839, -0.123158, 0.327244, -0.382768, 0.264523, -0.0261899, -0.224031, 0.372549, -0.351936, 0.171553, 0.0867126, -0.305612, 0.38577, -0.290796, 0.0638066, 0.19215, -0.360875, 0.365769, -0.204612, -0.0494312, 0.281037, -0.385057, 0.314269, -0.100809, -0.158417, 0.345724, -0.37608, 0.235703, 0.0116785, -0.253758, 0.380633, -0.334715, 0.13684, 0.123157, -0.327244, 0.382766, -0.264525, 0.0261936, 0.224035, -0.372546, 0.351939, -0.171553, -0.086716, 0.305609, -0.385773, ]
 ----------------------------------------------------------------
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|*                                                               |
|                                                                |
|                                                                |
|   *    *    *    *  *    *    *    *    *  *    *    *    *  * |
|     *          *      *               *      *          *      |
|          **                *     *                **           |
0 *    *        *        *    *   *    *        *        *    *  0
|  *                **                     **                *   |
|       *      *          *      *    *          *      *        |
|    *    *  *    *    *    *  *    *    *    *    *  *    *    *|
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
 ----------------------------------------------------------------


Impulse Response
[0.6, 0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ]
----------------------------------------------------------------
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|*                                                               |
|                                                                |
| *                                                              |
|                                                                |
|                                                                |
|                                                                |
0  **************************************************************0
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
|                                                                |
 ----------------------------------------------------------------
 Magnitude Response
 [1, 0.998844, 0.995378, 0.989612, 0.981561, 0.971248, 0.9587, 0.943952, 0.927044, 0.908025, 0.886946, 0.863869, 0.838861, 0.811995, 0.783354, 0.753026, 0.72111, 0.687715, 0.65296, 0.616979, 0.579924, 0.541968, 0.503315, 0.46421, 0.424957, 0.385947, 0.347699, 0.310931, 0.276655, 0.24631, 0.221863, 0.205697, 0.2, 0.205697, 0.221863, 0.24631, 0.276655, 0.310931, 0.347699, 0.385947, 0.424957, 0.46421, 0.503315, 0.541968, 0.579924, 0.616979, 0.65296, 0.687715, 0.72111, 0.753026, 0.783354, 0.811996, 0.838861, 0.863869, 0.886946, 0.908025, 0.927044, 0.943952, 0.9587, 0.971248, 0.981561, 0.989612, 0.995378, 0.998844, ]

-------------------------------NN-------------------------------
|*                                                               |
| *********                                             *********|
|          ****                                     ****         |
|              ***                               ***             |
|                 ***                         ***                |
|                    ***                   ***                   |
|                       **               **                      |
|                         ***         ***                        |
|                            *********                           |
|                                                                |
|                                                                |
 -------------------------------NN-------------------------------

 NN = Nyquist (half the sample rate)
*/
	#include <iostream>
	#include <complex>

	#define PI 3.141592653589793
	#define SAMPLE_RATE 64

	using namespace std;

	// This lowpass is a weighted moving average
	// If you change the + to - it will become a highpass filter
	void lowpass(float a, float b, float input, float output, int size){

	output[0] = a * input[0]; // avoid the negative index for n = 0

	for(int n = 1; n < size; n++){
	output[n] = a * input[n] + b * input[n - 1];
	}
	}


	// Print out the values of an array.
	void print_buffer(float *buffer, int size){
	cout << "[";
	for(int n = 0; n < size; n++){
	cout << buffer[n] << ", ";
	}
	cout << "]\n\n";
	}


	// Crappy O(n^2) discrete fourier transform code
	// ∑ e^jΩ * input[n]
	void dft(float input, std::complex<float> output, int size){
	float root = 2 * 3.141592653589793 * 1 / (float)size;
	for(int bin = 0; bin < size; bin++){
	complex<float> c = polar(1.f, root * bin);
	complex<float> sum = complex<float>(0, 0);
	for(int n = 0; n < size; n++){
	sum += input[n] * std::pow(c, n);
	}
	output[bin] = sum;
	}
	}


	// This is some kind of attempt to display a frequency response
	void print_frequency_response_graph(float *input, int size){
	const char *line = " -------------------------------NN-------------------------------";

	int quantized[size];
	for(int n = 0; n < size; n++){
	quantized[n] = (int)(input[n] * 10);
	}

	cout << line << endl;
	for(int magnitude = 10; magnitude >= 0; magnitude--){
	cout << "\|";
	for(int n = 0; n < size; n++){
	if(quantized[n] == magnitude){
	cout << "*";
	}else{
	cout << " ";
	}
	}
	cout << "\|" << endl;
	}
	cout << line << endl;
	}

	// This tries to print an audio buffer
	void print_audio_buffer(float *input, int size){
	const char *line = " ----------------------------------------------------------------";

	int quantized[size];
	for(int n = 0; n < size; n++){
	quantized[n] = (int)(input[n] * 10);
	}

	cout << line << endl;
	for(int magnitude = 10; magnitude >= -10; magnitude--){
	cout << (magnitude == 0 ? "0" : "\|");
	for(int n = 0; n < size; n++){
	if(quantized[n] == magnitude){
	cout << "*";
	}else{
	cout << " ";
	}
	}
	cout << (magnitude == 0 ? "0" : "\|") << endl;
	}
	cout << line << endl;
	}



	int main(void){

	float size = SAMPLE_RATE;

	// Change these:
	float frequency_hz = 25.f; // This is the frequency of the input cosine wave

	// Here are the coefficients
	// Make a + b = 1 so the filter has no overall gain
	float b = 0.4f;
	float a = 1.f - b;


	// Make a buffer holding 1 second of audio at the sample rate
	// and fill it up a 25hz sine wave
	float input[SAMPLE_RATE];
	float w = 2 * PI * (frequency_hz / (float)SAMPLE_RATE); // make this radian frequency

	for(int n = 0; n < size; n++){
	input[n] = cos(w * n);
	}

	// Buffer to hold the output of the lowpass filter
	float low[SAMPLE_RATE] = {0};

	// Buffer to hold an impulse, ie [1, 0, 0, 0, 0, 0, ...]
	float impulse[SAMPLE_RATE] = {0};
	impulse[0] = 1.f;

	// And a buffer to hold the filter's reponse to an impulse
	float impulse_response[SAMPLE_RATE] = {0};


	// Display the input buffer
	cout << "Input: " << endl;
	print_buffer(input, size);
	print_audio_buffer(input, size);

	// Low pass filter the input
	lowpass(a, b, input, low, size);

	// Display the lowpassed output
	cout << "Lowpass Output: " << endl;
	print_buffer(low, size);
	print_audio_buffer(low, size);

	// Figure out the impulse response of this filter
	// All we do is send an impulse into the filter and see what comes out
	lowpass(a, b, impulse, impulse_response, size);

	// Display the impulse response
	cout << "Impulse Response" << endl;
	print_buffer(impulse_response, size);
	print_audio_buffer(impulse_response, size);

	// A complex valued buffer to hold the filter's frequency response
	std::complex<float> frequency_response[SAMPLE_RATE];

	// Make a buffer to hold the magnitude reponse
	float magnitude_response[SAMPLE_RATE];

	// Magically, the fourier transform of the impulse reponse IS
	// the frequency response
	dft(impulse_response, frequency_response, size);

	// The frequency response contains the magnitude response if you
	// take the absolute value of each complex value
	// This is what we see as a graph on graphic EQs :)
	for(int n = 0; n < size; n++){
	magnitude_response[n] = abs(frequency_response[n]);
	}

	// Display the magnitude reponse. But remember that frequency
	// spectrums of real signals are mirror images about the nyquist (center)
	// frequency. Generally we only see the first half of this in programs
	// so a low pass filter is going to look like a dip down and then back up
	// in the mirrored half.
	cout << "Magnitude Response" << endl;
	print_buffer(magnitude_response, size);
	print_frequency_response_graph(magnitude_response, size);

	return 0;
	}

	/* Output

	Input:
	[1, -0.77301, 0.19509, 0.471397, -0.923879, 0.95694, -0.55557, -0.0980161, 0.707107, -0.995185, 0.83147, -0.290285, -0.382684, 0.881921, -0.980786, 0.634394, -7.35157e-07, -0.634393, 0.980785, -0.881922, 0.382685, 0.290284, -0.831469, 0.995185, -0.707106, 0.0980195, 0.555569, -0.956939, 0.923881, -0.4714, -0.195088, 0.773009, -1, 0.773011, -0.19509, -0.471397, 0.92388, -0.95694, 0.555574, 0.0980127, -0.707104, 0.995184, -0.831471, 0.290286, 0.382682, -0.881921, 0.980785, -0.634393, -1.60923e-06, 0.634389, -0.980784, 0.881923, -0.382687, -0.290282, 0.831464, -0.995185, 0.707113, -0.0980171, -0.555564, 0.956941, -0.923882, 0.471394, 0.195086, -0.773013, ]
	----------------------------------------------------------------
	\|* \|
	\| * * * * * * * * \|
	\| * * * * \|
	\| * * * * \|
	\| * * \|
	\| * * \|
	\| * * \|
	\| * * \|
	\| * * \|
	\| * * \|
	0 * * * * * * 0
	\| * * \|
	\| * * \|
	\| * * \|
	\| * * \|
	\| * * \|
	\| * * \|
	\| * * * *\|
	\| * * * * \|
	\| * * * * * * * * \|
	\| * \|
	----------------------------------------------------------------

	Lowpass Output:
	[0.6, -0.0638063, -0.19215, 0.360874, -0.365769, 0.204613, 0.0494342, -0.281038, 0.385057, -0.314268, 0.100808, 0.158417, -0.345724, 0.376079, -0.235703, -0.0116777, 0.253757, -0.380636, 0.334714, -0.136839, -0.123158, 0.327244, -0.382768, 0.264523, -0.0261899, -0.224031, 0.372549, -0.351936, 0.171553, 0.0867126, -0.305612, 0.38577, -0.290796, 0.0638066, 0.19215, -0.360875, 0.365769, -0.204612, -0.0494312, 0.281037, -0.385057, 0.314269, -0.100809, -0.158417, 0.345724, -0.37608, 0.235703, 0.0116785, -0.253758, 0.380633, -0.334715, 0.13684, 0.123157, -0.327244, 0.382766, -0.264525, 0.0261936, 0.224035, -0.372546, 0.351939, -0.171553, -0.086716, 0.305609, -0.385773, ]
	----------------------------------------------------------------
	\| \|
	\| \|
	\| \|
	\| \|
	\|* \|
	\| \|
	\| \|
	\| * * * * * * * * * * * * * * \|
	\| * * * * * * \|
	\| ** * * ** \|
	0 * * * * * * * * * * 0
	\| * * \|
	\| * * * * * * * \|
	\| * * * * * * * * * * * * * *\|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	----------------------------------------------------------------


	Impulse Response
	[0.6, 0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ]
	----------------------------------------------------------------
	\| \|
	\| \|
	\| \|
	\| \|
	\|* \|
	\| \|
	\| * \|
	\| \|
	\| \|
	\| \|
	0 **************************************************************0
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	\| \|
	----------------------------------------------------------------
	Magnitude Response
	[1, 0.998844, 0.995378, 0.989612, 0.981561, 0.971248, 0.9587, 0.943952, 0.927044, 0.908025, 0.886946, 0.863869, 0.838861, 0.811995, 0.783354, 0.753026, 0.72111, 0.687715, 0.65296, 0.616979, 0.579924, 0.541968, 0.503315, 0.46421, 0.424957, 0.385947, 0.347699, 0.310931, 0.276655, 0.24631, 0.221863, 0.205697, 0.2, 0.205697, 0.221863, 0.24631, 0.276655, 0.310931, 0.347699, 0.385947, 0.424957, 0.46421, 0.503315, 0.541968, 0.579924, 0.616979, 0.65296, 0.687715, 0.72111, 0.753026, 0.783354, 0.811996, 0.838861, 0.863869, 0.886946, 0.908025, 0.927044, 0.943952, 0.9587, 0.971248, 0.981561, 0.989612, 0.995378, 0.998844, ]

	-------------------------------NN-------------------------------
	\|* \|
	\| ******* *******\|
	\| ** ** \|
	\| * * \|
	\| * * \|
	\| * * \|
	\| \|
	\| * * \|
	\| ********* \|
	\| \|
	\| \|
	-------------------------------NN-------------------------------

	NN = Nyquist (half the sample rate)
	*/