Дискретное преобразование фурье

[Petrow]

Задача: есть косинусоида, например 18*cos(100*2*3.14*t+12) необходимо посчитать ее спектр, т.е. разложить в ряд Фурье и "увидеть" одну палку на значении 100(т.е. все остальные нули).....

вот программа....сначала массив(arrayA[i]) забивается значениями синусоиды...а затем считаю фурье по всем известной формуле дискретного преобразования....суть в том что выдает мне то же количество значений(300) но другие числа...хотя должно быть только одно! т.к. спектр косинуса-1 значение(палка)

#include "stdafx.h"
#include "math.h"
#include <iostream>
using namespace std;
float main(){
	////////////////////////////////////////////////////////////////////

	float  x,d,c,y;// х отвечает за действ. часть, у за мнимую. d и c введены для подсчета cos и sin(для удобства;)
	float  summa_D=0;//изначальная сумма для к-того элемента(палки), которая будет считать действительную часть
	float  summa_M=0;//изначальная сумма для к-того элемента(палки),которая будет считать мнимую часть
	float t=0;
	int  i,k; //номер палки(к), 
	float N=300; // всего отсчетов
	float arrayA[300];
	int n; // номер отсчета, в результате дискретизации
	setlocale(LC_ALL,"russian_russia");// 
	cout<<"Номер   | Действующие    |   Мнимые |"<<endl;// 
	for (i=0; i<3; i++){
        arrayA[i]=cos(100*2*3.14*t+12);
         	t=t+1;
	}
   for (k=0; k<N; k++)
   {
	   for (i=0; i<3;i++)
	   {
        d=cos((2*3,14*i*k)/N); //подсчитаем косинусную сост-ую(дейсвительную)
        c=sin((2*3,14*i*k)/N);//подсчитаем синусную сост-ую(т.е. мнимую)
		x=arrayA[i]*d;
	    y=-(c*arrayA[i]);
		summa_D=summa_D+x; // собсна суммируем что получилось=)))
        summa_M=summa_M+y;
	   }
   cout<<k<<"           "<<summa_D<<"             "<<summa_M<<"     "<<endl;
   }
return(0); //
}

[KobolD]

Ты помоему где то не там посморел формулу для дискретного преобразования, точнее не до конца дочитал, т.к. не среднегеометрического ни тангенса я у тебя не нашел.
Читай тут: http://slonpts.narod.ru/algo/math/furie/06.html

Вот пример на C# который я когда-то писал

public class FFT
    {
        const int N = 500; //количество разбиений
        double[] InnerArray = new double[N];
        Complex[] Furie = new Complex[N];

        public class Complex
        {
            public double Re;//Реальная часть
            public double Im;//Мнимая
            public double Amplitude;//Амплитуда для АЧХ
            public double Faza;//Фаза для ФЧХ
            public double Frecuensy;//Частота гармоники
        }

        private void Form1_Shown(object sender, EventArgs e)
        {
            for (int i = 0; i < N; i++)
                InnerArray[i] =5* Math.Sin(2 * Math.PI * i/100)*Math.Cos(2 * Math.PI * i / N);//задаем форму сигнала

            DPF(InnerArray);

        }

        void DPF(double[] Inner)
        {
            int N = Inner.Length;
            double Arg;

            for (int k = 0; k < N; k++)
            {
                Furie[k] = new Complex();
                for (int n = 0; n < Inner.Length; n++)
                {
                    Arg = 2 * Math.PI * k * n / N;
                    Furie[k].Re += Inner[n] * Math.Cos(Arg);
                    Furie[k].Im -= Inner[n] * Math.Sin(Arg);
                }
                Furie[k].Amplitude = (Math.Sqrt(Math.Pow(Furie[k].Re, 2) + Math.Pow(Furie[k].Im, 2))) / N;
                Furie[k].Faza = Math.Atan(Furie[k].Im / Furie[k].Re / Math.PI * 180);
                Furie[k].Frecuensy = ((N - 1) * (k));

            }
        }
    }

Собственно у тебя в масиве Furie в полях Amplitude-будет амплитуда i-той гармоники Faza-фаза и Frecuensy - частота гармоники.
В примерес синусойдой, в связи с ограниченой точностью вычисления не во всех остальных гармониках будет нулевая амплитуда, так что возможно надо будет дополнительно вводить фильтрацию, вроде того, что если амплитуда меньше 0,00001 то не учитываем эту гармонику.

[Petrow]

переделал....но вместо одного значения продолжают появляться 1000!!

#include "stdafx.h"
#include "math.h"
#include <iostream>
using namespace std;
float main(){
	////////////////////////////////////////////////////////////////////

	float  Re,Im,Arg,A,F,Fr;
	float  summa_Re=0;//изначальная сумма для к-того элемента(палки), которая будет считать действительную часть
	float  summa_Im=0;//изначальная сумма для к-того элемента(палки),которая будет считать мнимую часть
	float t=0;
	int  k; //номер палки(к), 
	float N=1000; // всего отсчетов
	float arrayA[1000];
	int n; // номер отсчета, в результате дискретизации
	setlocale(LC_ALL,"russian_russia");// 
	for (n=0; n<1000; n++){
        arrayA[n]=5*cos(2*3.14*800*t);
         	t=t+0.000010;
	}
   for (k=0; k<N; k++)
   {
	   for (n=0; n<1000;n++)
	   {
	    Arg=2*3.14*n*k/N;
        Re=cos(Arg)*arrayA[n];
		Im=sin(Arg)*arrayA[n];
		summa_Re=summa_Re+Re; 
        summa_Im=summa_Im+Im;
	   }
	   A=sqrt((summa_Re)*(summa_Re)+(summa_Im)*(summa_Im));
	   F=atan(summa_Im/summa_Re);
	   Fr=((N - 1)*(k));
   cout<<k<<"           "<<Fr<<"             "<<A<<"     "<<F<<"     "<<endl;
   }
return(0); 
}

[KobolD]

забыл амплитуду на количество отсчетов ноделить ( в твоем случае на 1000)
Я же говорю что "фильтрафию" надо делать

if (A>0.0000001)
       cout<<k<<"           "<<Fr<<"             "<<A<<"     "<<F<<"     "<<endl;

[Petrow]

ага) Спасибо !! получилось!!

[vadaa]

123456789test