Полный фурье анализ. Работает. Требуется дописать поиск максимумов с …

…заданной погрешностью и определение соотношения частот

Makefile

@@ -8,6 +8,7 @@ SRC += fft.c
 
 OBJ = $(SRC:.c=.o)
 
+
 all: $(TARGET)
 
 $(TARGET): obj
@@ -16,6 +17,7 @@ $(TARGET): obj
 	#cp $(TARGET) ../bin
 
 obj: $(SRC)
+	#$(CC) -Wall -std=c99 -c -g $(SRC) 
 	$(CC) -Wall -c -g $(SRC) 
 clean:
 	rm -f *.o $(TARGET) 

fft_an.c

@@ -1,99 +1,155 @@
 #include <stdio.h>
-#include <assert.h>
+#include <stdlib.h>
+#include <string.h> // for memcmp
+#include <stdint.h> // for int16_t and int32_t
 #include <math.h>
-#include <errno.h>
-
-// Структура, описывающая заголовок WAV файла.
-struct WAVHEADER
+#include "fft.h"
+
+struct wavfile
 {
-	// WAV-формат начинается с RIFF-заголовка:
-
-	// Содержит символы "RIFF" в ASCII кодировке
-	// (0x52494646 в big-endian представлении)
-	char chunkId[4];
-
-	// 36 + subchunk2Size, или более точно:
-	// 4 + (8 + subchunk1Size) + (8 + subchunk2Size)
-	// Это оставшийся размер цепочки, начиная с этой позиции.
-	// Иначе говоря, это размер файла - 8, то есть,
-	// исключены поля chunkId и chunkSize.
-	unsigned long chunkSize;
-
-	// Содержит символы "WAVE"
-	// (0x57415645 в big-endian представлении)
-	char format[4];
-
-	// Формат "WAVE" состоит из двух подцепочек: "fmt " и "data":
-	// Подцепочка "fmt " описывает формат звуковых данных:
-
-	// Содержит символы "fmt "
-	// (0x666d7420 в big-endian представлении)
-	char subchunk1Id[4];
-
-	// 16 для формата PCM.
-	// Это оставшийся размер подцепочки, начиная с этой позиции.
-	unsigned long subchunk1Size;
-
-	// Аудио формат, полный список можно получить здесь http://audiocoding.ru/wav_formats.txt
-	// Для PCM = 1 (то есть, Линейное квантование).
-	// Значения, отличающиеся от 1, обозначают некоторый формат сжатия.
-	unsigned short audioFormat;
-
-	// Количество каналов. Моно = 1, Стерео = 2 и т.д.
-	unsigned short numChannels;
-
-	// Частота дискретизации. 8000 Гц, 44100 Гц и т.д.
-	unsigned long sampleRate;
-
-	// sampleRate * numChannels * bitsPerSample/8
-	unsigned long byteRate;
-
-	// numChannels * bitsPerSample/8
-	// Количество байт для одного сэмпла, включая все каналы.
-	unsigned short blockAlign;
-
-	// Так называемая "глубиная" или точность звучания. 8 бит, 16 бит и т.д.
-	unsigned short bitsPerSample;
-
-	// Подцепочка "data" содержит аудио-данные и их размер.
-
-	// Содержит символы "data"
-	// (0x64617461 в big-endian представлении)
-	char subchunk2Id[4];
-
-	// numSamples * numChannels * bitsPerSample/8
-	// Количество байт в области данных.
-	unsigned long subchunk2Size;
-
-	// Далее следуют непосредственно Wav данные.
-};
-
-int main(int argc, char * argv[])
-{
-	FILE *file;
-	file = fopen("test.wav", "rb");
-	if (!file)
+	char    id[4];          // should always contain "RIFF"
+	int32_t totallength;    // total file length minus 8
+	char    wavefmt[8];     // should be "WAVEfmt "
+	int32_t format;         // 16 for PCM format
+	int16_t pcm;            // 1 for PCM format
+	int16_t channels;       // channels
+	int32_t frequency;      // sampling frequency
+	int32_t bytes_per_second;
+	int16_t bytes_by_capture;
+	int16_t bits_per_sample;
+	char    data[4];        // should always contain "data"
+	int32_t bytes_in_data;
+} __attribute__((__packed__));
+
+int is_big_endian(void) {
+	union {
+		uint32_t i;
+		char c[4];
+	} bint = {0x01000000};
+	return bint.c[0]==1;
+}
+
+int main(int argc, char *argv[]) {
+	char *filename=argv[1];
+	FILE *wav = fopen(filename,"rb");
+	struct wavfile header;
+
+	if ( wav == NULL ) {
+		fprintf(stderr,"Can't open input file %s\n", filename);
+		exit(1);
+	}
+
+	// read header
+	if ( fread(&header,sizeof(header),1,wav) < 1 ) {
+		fprintf(stderr,"Can't read input file header %s\n", filename);
+		exit(1);
+	}
+
+	// if wav file isn't the same endianness than the current environment
+	// we quit
+	if ( is_big_endian() ) {
+		if (   memcmp( header.id,"RIFX", 4) != 0 ) {
+			fprintf(stderr,"ERROR: %s is not a big endian wav file\n", filename); 
+			exit(1);
+		}
+	} else {
+		if (   memcmp( header.id,"RIFF", 4) != 0 ) {
+			fprintf(stderr,"ERROR: %s is not a little endian wav file\n", filename); 
+			exit(1);
+		}
+	}
+
+	if (   memcmp( header.wavefmt, "WAVEfmt ", 8) != 0 
+		|| memcmp( header.data, "data", 4) != 0 
+			) {
+		fprintf(stderr,"ERROR: Not wav format\n"); 
+		exit(1); 
+	}
+	if (header.format != 16) {
+		fprintf(stderr,"\nERROR: not 16 bit wav format.");
+		exit(1);
+	}
+	fprintf(stderr,"format: %d bits", header.format);
+	if (header.format == 16) {
+		fprintf(stderr,", PCM");
+	} else {
+		fprintf(stderr,", not PCM (%d)", header.format);
+	}
+	if (header.pcm == 1) {
+		fprintf(stderr, " uncompressed" );
+	} else {
+		fprintf(stderr, " compressed" );
+	}
+	fprintf(stderr,", channel %d", header.pcm);
+	fprintf(stderr,", freq %d", header.frequency );
+	fprintf(stderr,", %d bytes per sec", header.bytes_per_second );
+	fprintf(stderr,", %d bytes by capture", header.bytes_by_capture );
+	fprintf(stderr,", %d bits per sample", header.bytes_by_capture );
+	fprintf(stderr,", %d bytes in data", header.bytes_in_data );
+	fprintf(stderr,"\n" );
+
+	printf("\nchunk=%d\n",header.bytes_in_data/header.bytes_by_capture);
+	//показатель двойки массива
+	int p2=(int)(log2(header.bytes_in_data/header.bytes_by_capture));
+	printf("log2=%d\n",p2);
+	long int size_array=1<<(int)(log2(header.bytes_in_data/header.bytes_by_capture));
+	printf("size array=%ld \n", size_array);
+	if ( memcmp( header.data, "data", 4) != 0 ) { 
+		fprintf(stderr,"ERROR: Prrroblem?\n"); 
+		exit(1); 
+	}
+	fprintf(stderr,"wav format\n");
+
+
+
+//*****************************************************//
+	FILE * logfile;
+	logfile = fopen("test.txt", "w+");
+	if (!logfile)
 	{
-		printf("Failed open file, error");
+		printf("Failed open file, error\n");
 		return 0;
 	}
-
-	struct WAVHEADER header;
-
-	fread(&header, sizeof(struct WAVHEADER), 1, file);
-
-	// Выводим полученные данные
-	printf("Sample rate: %ld\n", header.sampleRate);
-	printf("Channels: %d\n", header.numChannels);
-	printf("Bits per sample: %d\n", header.bitsPerSample);
-
-	// Посчитаем длительность воспроизведения в секундах
-	float fDurationSeconds = 1.f * header.subchunk2Size / (header.bitsPerSample / 8) / header.numChannels / header.sampleRate;
-	int iDurationMinutes = (int)floor(fDurationSeconds) / 60;
-	fDurationSeconds = fDurationSeconds - (iDurationMinutes * 60);
-	printf("Duration: %02d:%02.f\n", iDurationMinutes, fDurationSeconds);
-
-	fclose(file);
-
-	return 0;
+//*****************************************************//	
+
+	//fft!!!
+	float *c;  // массив поворотных множителей БПФ
+	float *in; // входной массив
+	float *out;// выходной массив
+
+	c = calloc(size_array*2, sizeof(float));
+	in = calloc(size_array*2, sizeof(float));
+	out = calloc(size_array*2, sizeof(float));
+
+	fft_make(p2,c);// функция расчёта поворотных множителей для БПФ
+
+	int16_t value;
+	int i=0;
+	int j=0;
+
+	while( fread(&value,sizeof(value),1,wav) ) {
+
+		//fprintf(logfile,"%d %hi\n",i, value);
+		in[j]=(float)value;
+		i++;
+		j+=2;
+		if (i > size_array)  break;
+	}
+	//fft_calc(p2, c,	in,	out, 0);
+	fft_calc(p2, c,	in,	out, 1);
+	j=0;
+	double delta=((float)header.frequency)/(float)size_array;
+	//printf("delta=%3.10f",delta);
+	double cur_freq=0;
+	for(i=0;i<(2*size_array);i+=2) {
+		//fprintf(logfile,"%.6f %f\n",cur_freq, (sqrt(out[i]*out[i]+out[i+1]*out[i+1])));
+		fprintf(logfile,"%.6f %f\n",cur_freq, (sqrt(out[i]*out[i]+out[i+1]*out[i+1])));
+		cur_freq+=delta;
+		//j++;
+	}
+
+	fclose(logfile);
+	fclose(wav);
+	//printf("%lld\n",sum);
+	exit(0);
 }

gnuplot.graph

@@ -5,12 +5,15 @@ set output "result.ps"
 #set output "result.png"
 #set multiplot layout 2,1
 set grid xtics ytics
-#set xlabel "N" 
-#set ylabel "Counter"
-set xrange [0:420]
+set log x
+set xlabel "Freq, Hz" 
+set ylabel "Amp"
+set xrange [10:22050]
 #set yrange [-3:3]
 #set style line 1 lt 1 lw 50
 #set style line 1 lt 1 lw 1
+plot "test.txt" using 1:2 title "AFC" with lines linestyle 1
+
 #plot "sounddist.dat" using 1:2 title "Распределение генератора случайных чисел" with boxes fs pattern 1, "sounddist.dat" using 1:3 title "Keybd stats" with lines linestyle 3
 
 #plot "result.dat" using 1:2 title "Distribution" with boxes fs pattern 1, "result.dat" using 1:3 title "Teor" with lines linestyle 3
@@ -21,7 +24,7 @@ set xrange [0:420]
 #plot "disp.dat" using 1:2 title "Disp" with boxes fs pattern 1
 
 #set style line 1 lt 1 pt 0
-plot "test.txt" using 1:2 with lines linestyle 1
+
 #set xrange [0:100]
 #plot 1543431*3.66*0.000001*exp(-3.66*0.000001*x)
 #plot 1543431*0.16666666*exp(-0.1666666*x)