Skip to content

Instantly share code, notes, and snippets.

@safiire
Created July 8, 2014 04:09
Show Gist options
  • Save safiire/03ce7a3e719617676176 to your computer and use it in GitHub Desktop.
Save safiire/03ce7a3e719617676176 to your computer and use it in GitHub Desktop.
One Zero Low Pass Filter
#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)
*/
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment