參與評(píng)論

評(píng)論

登錄后參與評(píng)論

全部評(píng)論(1條)

• 

  我昱云南的人民 2016-03-30 03:47:15

  這個(gè)可比你想象的復(fù)雜多了，s是個(gè)復(fù)變量，1/(s+1)極點(diǎn)在-1，要想用C語(yǔ)言寫(xiě)，必須理解清楚下面幾個(gè)問(wèn)題： 1、輸入必須是個(gè)有限序列，比如(x+yi)，x和y分別是兩個(gè)長(zhǎng)度為N的數(shù)組 2、要過(guò)濾的頻率，必須是個(gè)整型值，或者是個(gè)整型區(qū)間 3、輸出結(jié)果同樣是兩個(gè)長(zhǎng)度為N的數(shù)組(p+qi) 4、整個(gè)程序需要使用Z基本的復(fù)數(shù)運(yùn)算，這一點(diǎn)C語(yǔ)言本身不提供，必須手工寫(xiě)復(fù)函數(shù)運(yùn)算庫(kù) 5、實(shí)現(xiàn)的時(shí)候具體算法還需要編，這里才是你問(wèn)題的核心。 我可以送你一段FFT的程序，自己琢磨吧，和MATLAB的概念差別很大： #include <assert.h> #include <math.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <windows.h> #include "complex.h" extern "C" { // Discrete Fourier Transform (Basic Version， Without Any Enhancement) // return - Without Special Meaning， constantly， zero int DFT (long count， CComplex * input， CComplex * output) { assert(count); assert(input); assert(output); CComplex F， X， T， W; int n， i; long N = abs(count); long Inversing = count < 0? 1: -1; for(n = 0; n < N ; n++){ // compute from line 0 to N-1 F = CComplex(0.0f， 0.0f); // clear a line for(i = 0; i < N; i++) { T = input[i]; W = HarmonicPI2(Inversing * n * i， N); X = T * W; F += X; // fininshing a line }//next i // save data to outpus memcpy(output + n， &F， sizeof(F)); }//next n return 0; }//end DFT int fft (long count， CComplex * input， CComplex * output) { assert(count); assert(input); assert(output); int N = abs(count); long Inversing = count < 0? -1: 1; if (N % 2 || N < 5) return DFT(count， input， output); long N2 = N / 2; CComplex * iEven = new CComplex[N2]; memset(iEven， 0， sizeof(CComplex) * N2); CComplex * oEven = new CComplex[N2]; memset(oEven， 0， sizeof(CComplex) * N2); CComplex * iOdd = new CComplex[N2]; memset(iOdd ， 0， sizeof(CComplex) * N2); CComplex * oOdd = new CComplex[N2]; memset(oOdd ， 0， sizeof(CComplex) * N2); int i = 0; CComplex W; for(i = 0; i < N2; i++) { iEven[i] = input[i * 2]; iOdd [i] = input[i * 2 + 1]; }//next i fft(N2 * Inversing， iEven， oEven); fft(N2 * Inversing， iOdd， oOdd ); for(i = 0; i < N2; i++) { W = HarmonicPI2(Inversing * (- i)， N); output[i] = oEven[i] + W * oOdd[i]; output[i + N2] = oEven[i] - W * oOdd[i]; }//next i return 0; }//end FFT void __stdcall FFT( long N， // Serial Length， N > 0 for DFT， N < 0 for iDFT - inversed Discrete Fourier Transform double * inputReal， double * inputImaginary， // inputs double * AmplitudeFrequences， double * PhaseFrequences) // outputs { if (N == 0) return; if (!inputReal && !inputImaginary) return; short n = abs(N); CComplex * input = new CComplex[n]; memset(input， 0， sizeof(CComplex) * n); CComplex * output= new CComplex[n]; memset(output，0， sizeof(CComplex) * n); double rl = 0.0f， im = 0.0f; int i = 0; for (i = 0; i < n; i++) { rl = 0.0f; im = 0.0f; if (inputReal) rl = inputReal[i]; if (inputImaginary) im = inputImaginary[i]; input[i] = CComplex(rl， im); }//next i int f = fft(N， input， output); double factor = n; //factor = sqrt(factor); if (N > 0) factor = 1.0f; else factor = 1.0f / factor; //end if for (i = 0; i < n; i++) { if (AmplitudeFrequences) AmplitudeFrequences[i] = output[i].getReal() * factor; if (PhaseFrequences) PhaseFrequences[i] = output[i].getImaginary() * factor; }//next i delete [] output; delete [] input; return ; }//end FFT int __cdecl main(int argc， char * argv[]) { fprintf(stderr， "%s usage:\n"， argv[0]); fprintf(stderr， "Public Declare Sub FFT Lib \"wfft.exe\" \ (ByVal N As Long， ByRef inputReal As Double， ByRef inputImaginary As Double， \ ByRef freqAmplitude As Double， ByRef freqPhase As Double)"); return 0; }//end main };//end extern "C"

  贊(15)

  回復(fù)(0)

  評(píng)論

  評(píng)論

  登錄后參與評(píng)論