Trouble Generating a Sine Wave in C/Matlab from Basics

Question 1

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;
}

Question 2

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.

Question 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;
}