Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Алгоритмы > точность Фурье

Автор: stab 1.2.2003, 20:54

Сгенерял я сигнал с частотой 1000Гц, сделал Фурье и получил и что в сигнале присутствуют много разных частот, от 500 и до 1500 (примерно), причем амплитуды у них довольно большие. Делаю опять Фурье получаю вроде бы исходный сигнал, но он совершенно не похож на оригинал. Подскажите как быть?

Автор: podval 2.2.2003, 05:55

А какой именно сигнал и как ты сгенерил?

Автор: stab 3.2.2003, 21:26

синусойда, сгенерял в CoolEdit (редактор wav-ников если кто незнает)

Автор: podval 4.2.2003, 01:45

Чистая синусоида - это не сигнал. Она не несет информации. Сигнал же несет информацию. А информация - это какие-либо нерегулярности. Допустим, модулированная синусоида - это уже сигнал. Чем-нибудь модулировал синусоиду? Чистая она или нет? Какой алгоритм преобразования применял? БПФ или в лоб ДПФ?
Такой уж я буквоед

Автор: stab 4.2.2003, 04:59

немодулированая, просто синусоида, ДПФ. Про сигнал: ну лана не сигнал, а отсчеты

Автор: podval 4.2.2003, 07:06

Взглянуть бы своими глазами. А так все равно, что вырезать аппендицит по телефону.
Ну я понимаю, если бы БПФ с каким-то окном. Там имеется "натяжка": предполагается, что за пределами окна сигнал продолжается периодически. Отсюда некоторый гемморой.
Ну а у тебя что? Последствия дискретизации, квантования?
Можешь как-нибудь нарисовать или кинуть код?

Автор: val 20.2.2003, 23:46

Вот пример расчетая прямого и обратного преобразования Фурье... Проверь, как у тебя...
#define Pi 3.141592
#define Window 160
#define HalfWin 80
double rex[HalfWin+1];
double imx[HalfWin+1];
double workbuf[Window]; //тут лежит твой сигнал

void DFT()//Discret Fourier Transform function
{
unsigned int i;
unsigned int k;
double w;
w=2*Pi/Window;
for(k=0;k<=HalfWin;k++)
{
rex[k]=0;
imx[k]=0;
}
for(k=0;k<=HalfWin;k++)
for(i=0;i<Window;i++)
{
rex[k]=rex[k]+workbuf[i]*cos(w*k*i);
imx[k]=imx[k]-workbuf[i]*sin(w*k*i);
}

void IDFT()//Inverse Discret Fourier Transform function
{
for(k=0;k<=HalfWin;k++)
{
rex[k]=rex[k]/HalfWin;
imx[k]=-imx[k]/HalfWin;
}
rex[0]=rex[0]/2;
rex[HalfWin]=rex[HalfWin]/2;

for(i=0;i<Window;i++)
workbuf[i]=0;

for(k=0;k<=HalfWin;k++)
for(i=0;i<Window;i++)
{
workbuf[i]=workbuf[i]+rex[k]*cos(w*k*i);
workbuf[i]=workbuf[i]+imx[k]*sin(w*k*i);
}
}
void main()
{
Ввести массив workbuf;
DFT();
IDFT();
//рассчитали прямое и обратное Фурье… Сигнал не изменился
}

Автор: podval 21.2.2003, 04:34

На, поиграйся с моим БПФ

Код


/*
Основная функция, реализующая БПФ.
sig - массив входных данных, результат возвращается сюда же.
n - размер массива sig. Должен быть степенью двойки. Недостающие
    отсчёты дополняются нулями.
sn - массив размером n+1, представляющий ядро преобразования.
p = Log2(n) или n = pow(2,p), поэтому в принципе необязательно
    обе величины задавать в аргументах функции, это надо было
    в учебных целях и еще это в вычислительном смысле удобно.
*/
//---------------------------------------------------------------------------
void fft(complex<double> *sig, int p, complex<double> *sn, int n)
{
  complex<double> *spr = new complex<double>[n];
  int r2 = 1;
      r2 <<= (p - 1);
  int l1 = r2 + r2;
  for(int ir=1; ir <= p; ir++)
    {
     int i1 = p - ir + 1;
     int b = 1;
  b <<= (i1 - 1);
     int r = 1;
  r <<= (p - i1);
  r -= 1;
     int rm = r + 1;
     for(int mr = 1; mr <= rm; mr++)
       {
	int m = mr - 1;
	int r1 = b*m;
	int r3 = b*(m + r + 1);
	int l = r1 + r1;
	int m31 = r1 + 1;
	int m41 = r3 + 1;
	for(int km = 1; km <= b; km++)
   {
    int k = km - 1;
    int m11 = k + l + 1;
    int m21 = m11 + b;
    int m51 = k + r1 + 1;
    int m61 = k + r2 + r1 + 1;
    spr[m51-1] = sig[m11-1] + sig[m21-1]*sn[m31-1];
    spr[m61-1] = sig[m11-1] + sig[m21-1]*sn[m41-1];
   }
       }
     for(int m = 0; m < l1; m++)
         sig[m] = spr[m];
    }
   for (int i=0; i < n;i++)
        sig[i] /= n;
delete[] spr; spr = 0;
}
//--------------------------------------------------------------------
/*
Пример работы с функцией fft
Здесь предполагается, что имеются некоторые данные, а
юзер хочет сделать БПФ применительно к выборке между отсчётами
с номерами от FFTBegin до FFTEnd
*/
   // в следующем цикле do...while готовятся данные для БПФ
   // а именно размер входного массива и показатель степени двойки

   int p1, n = 2, nn = 1;   // p1 = 2 в степени nn
   do {
      p1 = (n <<= 1); // p1 - сюда пишем размер преобразуемого массива
      nn += 1;        // nn - показатель степени двойки (p1 = pow(2,nn))
   } while((FFTEnd - FFTBegin) > p1);

   complex<double> *sig = new complex<double>[p1]; // массив входных данных
   complex<double> *vyx = new complex<double>[p1+1]; // ядро преобразования

   // готовим ядро преобразования
   double a = 2*M_PI/p1;
   double b1;
   for(int i=0; i < (p1+1); i++)
   {
      b1 = a*i;
      vyx[i] = complex<double>(cos(b1),sin(b1));
   }

   // здесь вставить ввод исходных данных в массив sig
   // обращаю внимание, что он complex!!!
   // это значит, что если данные представлены вещественными
   // числами, их надо привести к комплексному виду с нулевой
   // мнимой частью.

   // процедура БПФ
   fft(sig, nn, vyx, p1);


   // можно сделать и обратное преобразование, для этого надо слегка 
   // переделать ядро (найти комплексно-сопряженные числа)

   for(int i=0; i < (p1+1); i++)
       vyx[i] = conj(vyx[i]);

   // а затем вызвать fft

/*
NB: для работы необходимо
#include<complex.h>
#include<math.h>
*/

Автор: Crait 21.2.2003, 09:32

А выборку на оконную функцию домножать не пробовали ? -
есть такие, имеющие колоколообразную форму.
Используйте, скажем, окно Хэмминга - из наиболее известных.
f = 0.54 - 0.46 * cos(x)

Автор: podval 21.2.2003, 16:47

Цитата

А выборку на оконную функцию домножать не пробовали ?

И не только на Хэмминга

Автор: Iskanderus 1.3.2007, 07:29

Уважаемый podval, в ваши функции исходный сигнал передается в виде
complex<double> *sig

Обяъсните, пожалуйста, каким образом этот сигнал превращается в массив комплексных чисел типа double?

Автор: podval 2.3.2007, 14:54

В функцию передали указатель на массив. Далее делаем, что хотим.