Trouble Generating a Sine Wave in C/Matlab from Basics

Question

I've been writing a basic C program to generate a sine wave and write it out to STDOUT for piping into a filter which I am also going to write but I'm having issues.

I am able to generate what sounds like a good clean sine wave, but the frequency is off. When I listen to the output I am generating and compare to a true signal generators output, the frequency of mine is slightly too high.

The C code should generate raw 16 bit signed stereo wav @ 44100Hz for a frequency passed as an argument to the program, or 500Hz as default.

A matlab version of the code is here, (modded very slightly for 1 indexing in Matlab, WHY MATHWORKS WHY!?) minus the passing to STDOUT etc as I know that works fine

CHANNELS = 2; SAMPLING_RATE = 44100; NUM_SAMPLES = 512;

frequency = 1000;

%% BEGIN WHILE(1) LOOP HERE

output_buff = zeros(1,CHANNELS*NUM_SAMPLES);

for i = 1:2:CHANNELS*NUM_SAMPLES
    output_buff(i) = 30000 * sin(frequency * pi * (i-1)/NUM_SAMPLES);
    output_buff(i+1) = output_buff(i);
end

%% OUTPUT TO STDOUT AND REPEAT

I should add that this code runs (in the C version) inside a while true loop, generating a full output_buff of values, then pushing the buffer to STDOUT.

I have written some further testing code to see what's actually being generated as follows:

plot(1:CHANNELS*NUM_SAMPLES, output_buff)

output_buff = output_buff .* hanning(length(output_buff))';    
Y = fft(output_buff);    
Mag=abs(Y(1:length(output_buff)/2)).^2;    
[a,b]=max(Mag);    
% Result    
SAMPLING_RATE*b/length(output_buff)

When I run this script, I can see at the end the frequency of the generated signal is actually 1.0767e+03Hz... Close but no cigar...

I've tried tweaking some of the parameters but I have no idea what's wrong or how to make the generated frequency more accurate.

The C code itself is on my Linux install and I can add that tomorrow if needs be.

Answer 1

You are dealing with Digital Signal Processing, it's a complex field, and there is a dsp stackexchange site dedicated to it. My advises are:

• If you want to output something which is exactly what you entered, choose to generate a frequency that is a division of the sample rate by an power of two, for instance 44100/44 = 1002.272727... In that case, a sample that will be a power of two long, will perfectly fit in your FFT input.

• If you want better FFT results, try increasing the number of samples to 4096 or 8192. Because 512 samples at a sample rate of 44,100Hz means you have 512/44,100=~ 11.61 ms of signal..so it means you have an incomplete numbers of sine waves. One complete sine wave at 1000 Hz is exactly 1 ms. This incomplete number of periods could lead to approximation errors by the FFT.

Answer 2

So turns out I had a few things wrong. My for loop was using the number of bytes in the output buffer, not the number of elements in the buffer. I was overlooking the relevance of the wav file sampling rate in generating the sine wave, and the discontinuity in using the loop variable was causing strange artifacts.

The final working code is as follows:

#include <sys/types.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>

#include <unistd.h>

#define M_PI 3.1415926535897932384626433832795L
#define NUM_CHANNELS 2
#define NUM_SAMPLES 512
#define SAMPLE_RATE 44100

double frequency = 500.0;

int main(int argc, char *argv[])
{   
    if (argc >= 2){
        frequency = atof(argv[1]);
    }

    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);
    int counter = 0;

    while(1){
        for (int i = 0; i < NUM_CHANNELS * NUM_SAMPLES; i+=NUM_CHANNELS){
            samples[i] = (int16_t) 3000 * sin(frequency * 2 * M_PI * (double)counter /(double)(SAMPLE_RATE));
            samples[i+1] = samples[i];
            counter++;
        }

        int done=write(STDOUT_FILENO, samples, cbBuffer);
        if(done<0){
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }else if(done!=cbBuffer){
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}

Answer 3

I think I have a slightly improved version - using a quarter sine-lookup table which is pre-defined outside the main while-loop.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <stdint.h>
#include <unistd.h>

// PI defined here for use in code
#define PI 3.141592653589793
// SINE_TABLE_SIZE - number of values in the sine lookup table - defined here for use in code
#define SINE_TABLE_SIZE 256
#define NUM_CHANNELS 2
#define NUM_SAMPLES 512

//* Sampling frequency in HZ. Must only be set to 8000, 16000, 24000 32000, 44100 (CD standard), 48000 or 96000 */
int sampling_freq = 44100;
// Array of data used by sinegen to generate sine. These are the initial values.
float y[3] = {0,0,0};
float x[1] = {1}; // impulse to start filter
float a0 = 1.4142; // coefficients for difference equation
float b0 = 0.707;
// Holds the value of the current sample
float sample;
float sine_freq = 500.0;

// An array of floats containing SINE_TABLE_SIZE elements
float table[SINE_TABLE_SIZE];
// Step variable for use in the quarter sinegen function
float step4 = 0;
// Modified step value to within Sine wave table values - for use in quarter sine wave table
float table_value;



/********************************** sine_init() **************************/ 
void sine_init4()
{
    // Fill the lookup table with 256 sine data points across one quarter cycle.
    int j;
    for(j=0; j < SINE_TABLE_SIZE; j++)
    {
        table[j] = sin(((0.5*PI*j)/SINE_TABLE_SIZE));
    }
}



/********************************** sinegen4() **************************/
float sinegen4(void)
{
/* This code produces a variable frequency sine wave, using a
quarter-sine-wave lookup table.
*
* */
    float wave4 = 0; // Initialise a variable to store the sine wave datapoint values in.
// To create a sine wave from the quarter sinewave table data.
//For values in the first sinewave quadrant - no adjustment to the step value needs to be made.
    if (step4 < (SINE_TABLE_SIZE))
    {
        table_value = step4;
        wave4 = table[(int)step4];
    }
//Second quadrant - step value must be adjusted to bring the value back into the range 0-255
    else if (step4 < (2*SINE_TABLE_SIZE) && (step4 >= SINE_TABLE_SIZE))
    {
        table_value = ((SINE_TABLE_SIZE-1)-(step4-SINE_TABLE_SIZE));

        wave4 = table[(int)((SINE_TABLE_SIZE-1)-(step4- SINE_TABLE_SIZE))];
    }
//Third quadrant - step value must be adjusted to bring the value back into the range 0-255 and the wave value negated
    else if (step4 < (3*SINE_TABLE_SIZE) && (step4 >= (2*SINE_TABLE_SIZE)) )
    {   
        table_value = (step4-(2*SINE_TABLE_SIZE));
        wave4 = -table[(int)(step4-(2*SINE_TABLE_SIZE))];
    }
//Fourth quadrant - step value must be adjusted to bring the value back into the range 0-255 and the wave value negated
    else if (step4 < (4*SINE_TABLE_SIZE) && (step4 >=(3*SINE_TABLE_SIZE)) )
    {
        table_value = ((SINE_TABLE_SIZE-1)-(step4-(3*SINE_TABLE_SIZE)));
        wave4 = -table[(int)((SINE_TABLE_SIZE-1)-(step4-(3*SINE_TABLE_SIZE)))];
    }
// assign step a value based on sampling frequency and desired output frequency to calculate next table value required.
    step4 += ((4*SINE_TABLE_SIZE)/(sampling_freq/sine_freq));
//To prevent step containing values greater than 4*SINE_TABLE_SIZE-1 which would cause the operation to overflow.
    if (step4 > ((4*SINE_TABLE_SIZE-1)))
    {
        step4 = step4 - (4*SINE_TABLE_SIZE-1);
    }

    return wave4;
}




int main(int argc, char *argv[])
{   

    if(argc > 1)
    {
        sine_freq = atoi(argv[1]);
//      printf("n = %d \n", n );
    }

    // initialises table of one quarter sinewave data
    sine_init4();
    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);
    // Loop endlessley generating a sine wave
    while(1)
    {

    // Calculate next sample from quarter sinewave data table
        for (int i = 0; i < NUM_CHANNELS * NUM_SAMPLES; i+=NUM_CHANNELS)
        {
            samples[i] = 3000*sinegen4();
            samples[i+1] = samples[i];
            //printf(" samples[i] = %d", samples[i]);
        }
        // Copy one sample to output
        int done=write(STDOUT_FILENO, samples, cbBuffer);
        //int done = cbBuffer;
        if(done<0)
        {
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }
        else if(done!=cbBuffer)
        {
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}