从基础知识在 C/Matlab 中生成正弦波的麻烦答案

【问题标题】：Trouble Generating a Sine Wave in C/Matlab from Basics从基础知识在 C/Matlab 中生成正弦波的麻烦
【发布时间】：2014-02-04 22:47:22
【问题描述】：

我一直在编写一个基本的 C 程序来生成一个正弦波并将其写入 STDOUT 以管道到一个过滤器中，我也将编写该过滤器，但我遇到了问题。

我能够产生听起来很干净的正弦波，但频率已关闭。当我听到我正在生成的输出并与真正的信号发生器输出进行比较时，我的频率有点太高了。

C 代码应为作为参数传递给程序的频率或默认为 500Hz 的频率生成原始 16 位有符号立体声 wav @ 44100Hz。

代码的 matlab 版本在这里，（对 Matlab 中的 1 个索引进行了非常轻微的修改，WHY MATHWORKS WHY！？）减去传递给 STDOUT 等，因为我知道它工作正常

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

我应该补充一点，这段代码在 while true 循环中运行（在 C 版本中），生成一个完整的 output_buff 值，然后将缓冲区推送到 STDOUT。

我已经编写了一些进一步的测试代码，以查看实际生成的内容如下：

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)

当我运行这个脚本时，我可以看到最后生成信号的频率实际上是 1.0767e+03Hz... Close but no cigar...

我已经尝试调整一些参数，但我不知道出了什么问题或如何使生成的频率更准确。

C 代码本身在我的 Linux 安装中，如果需要，我可以在明天添加。

【问题讨论】：

你的音频输出设备是什么？你知道它的时钟是准确的吗？附：因为 FORTRAN
在我的联想 Y510p、3.2GHz i7、16GB RAM 上运行代码，通过系统默认音频系统运行

标签： c matlab audio wav trigonometry

【解决方案1】：

您正在处理数字信号处理，这是一个复杂的领域，并且有一个专门的 dsp stackexchange 站点。我的建议是：

如果您想输出与您输入的内容完全相同的内容，请选择生成一个频率，该频率是采样率除以 2 的幂，例如 44100/44 = 1002.272727... 在这种情况下，一个长为 2 的幂的样本将完全适合您的 FFT 输入。
如果您想要更好的 FFT 结果，请尝试将样本数增加到 4096 或 8192。因为 44,100Hz 的采样率下的 512 个样本意味着您有 512/44,100=~ 11.61毫秒的信号..所以这意味着您的正弦波数量不完整。一个完整的 1000 Hz 正弦波正好是 1 ms。这种不完整的周期数可能会导致 FFT 的近似误差。

【讨论】：

【解决方案2】：

原来我做错了一些事情。我的 for 循环使用的是输出缓冲区中的字节数，而不是缓冲区中的元素数。我忽略了 wav 文件采样率与生成正弦波的相关性，而使用循环变量的不连续性导致了奇怪的伪影。

最终的工作代码如下：

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

【讨论】：

【解决方案3】：

我想我有一个稍微改进的版本 - 使用在主 while 循环之外预定义的四分之一正弦查找表。

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

【讨论】：