http://ru.wikipedia.org/wiki/%D0%9E%D0%BA%D0%BE%D0%BD%D0%BD%D0%BE%D0%B5_%D0%BF%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A4%D1%83%D1%80%D1%8C%D0%B5#.D0.9E.D0.BA.D0.BD.D0.BE_.D0.A5.D0.B0.D0.BD.D0.BD.D0.B0 http://ru.wikipedia.org/wiki/%D0%94%D0%B8%D1%81%D0%BA%D1%80%D0%B5%D1%82%D0%BD%D0%BE%D0%B5_%D0%BF%D1%80%…

gistfile1.cpp

#define _CRT_SECURE_NO_DEPRECATE
#include <stdio.h>
#include <conio.h>
#include <iostream>
#include <math.h>
#include <time.h>
#include <fftw3.h>

#define M_PI 3.1415926535897932384626433832795


int M = 4;
const int dimArr = int(pow(2.0, M));	//Размер массива, оно же количество частей, на которое мы разбиваем сигнал. Должно быть степенью 2ки.
const double w = 4.0; // Частота w

using namespace std;

double function(double x)	//Заданный сигнал
{
	return sin(w*x);
}
double rectangle (double y) //Прямоугольное окно
{
	if(y>=0 && y <= (double)dimArr)
		return 1;
	else 
		return 0;
}

double hanning_window (double x)	//Окно Ханна
{
	if (x <= (double)dimArr) 	
		return 0.5*(1-cos(2.0*M_PI*x/double((double(dimArr)-1))));
	else return 0;
}

void main (void) {
	int i;
	double dx = 1.0/dimArr;
	double xCur;
	double x=0;

	FILE *file = fopen("file.txt", "w");	//Сюда запишется спектральная мощность
	FILE *file2 = fopen("file2.txt", "w");	//Сюда запишем сигнал, после обратного преобразования Фурье

	double * fx0, *fx1;	//Массивы для действительной (0) и мнимой (1) части входного сигнала
	double * fk0, *fk1;	//Тоже самое для Спектра входного сигнала


	fx0 = (double*)malloc(sizeof(double) * dimArr);
	fx1 = (double*)malloc(sizeof(double) * dimArr);

	fk0 = (double*)malloc(sizeof(double) * dimArr);
	fk1 = (double*)malloc(sizeof(double) * dimArr);
	
	for (i=0; i<dimArr; i++)		//Задание начальных условий
	{
		xCur = 2.0*M_PI*i*dx;
		fx0[i] = function(xCur)*hanning_window((double)i);	//Сразу промоделируем входной сигнал - в данном случае окном Ханна
		//fx0[i] = function(xCur)*rectangle(i);	//А в этом - прямоугольником
		fx1[i] = 0.0;
	}

	for (int k = 0; k < dimArr; k++)	//Делаем прямое преобразование Фурье
	{
		fk0[k] = 0.0;
		fk1[k] = 0.0;
		for (i = 0; i < dimArr; i++)
		{
			fk0[k] += fx0[i]*cos(2*M_PI*i*k/dimArr) - fx1[i]*sin(2*M_PI*i*k/dimArr);
			fk1[k] += fx0[i]*sin(2*M_PI*i*k/dimArr) + fx1[i]*cos(2*M_PI*i*k/dimArr);
		}
	}

	for (i=0; i<dimArr; i++)	//Необходимо также пронормировать спектр
	{
		fk0[i] = fk0[i]/(double)dimArr;
		fk1[i] = fk1[i]/(double)dimArr;
	}

	for (i=0; i<dimArr; i++)	
	{
		fk0[i] = (fk0[i]*fk0[i] + fk1[i]*fk1[i])/(2*M_PI);	//По задаче нам интересна Спектральная Мощность - т.е. полная амплитуда. Это Re*Re + Im*Im от комплексного значения компоненты спектра
		fprintf(file, "%d %1.5f\n", i, fk0[i]);		//Запись спектра в Файл
	}

	for (int k = 0; k < dimArr; k++)	//Обратное преобразование Фурье (должно получить исходный сигнал)
	{
		fx0[k] = 0.0;
		fx1[k] = 0.0;
		for (i = 0; i < dimArr; i++)
		{
			fx0[k] += fk0[i]*cos(2*M_PI*i*k/dimArr) + fk1[i]*sin(2*M_PI*i*k/dimArr);
			fx1[k] += -fk0[i]*sin(2*M_PI*i*k/dimArr) + fk1[i]*cos(2*M_PI*i*k/dimArr);
		}
	}

	
	for (i=0; i<dimArr; i++)	//Запись сигнала после всех процедур в Файл
	{
		fprintf(file2, "%1.5f %1.5f\n", 2*M_PI*i/dimArr, fx0[i]);
	}

	delete(fx0);
	delete(fx1);
	delete(fk0);
	delete(fk1);

	fclose(file);
	fclose(file2);

}