I am trying to write a code in c++ equivalents to the following Matlab code:
n = 32;
s1 = dlmread('ssa.txt');
s2 = dlmread('ssb.txt');
fd1 = fft(s1);
fd2 = fft(s2);
nf = sqrt((fd1 * fd1') * (fd2 * fd2')) / n
Where
ssa.txt= 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
and
ssb.txt= 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 1
Unfortunately the result I am getting in c++ is bit different from Matlab
result from Matlab:
nf = 12
result in c++:
nf(4.88621,0)
My c++ code :
//iterDelayEst.cpp
#include <complex>
#include "binFreq.cpp"
#include <valarray>
typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;
double iterDelayEst(int n, CArray& x, CArray& y)
{
/**********************************************constants************************************************
*******************************************************************************************************/
//exit if uncertainty below threshold
fft(x);
fft(y);
auto fd2Tau = y;
//frequency domain representation of signals
std::vector<double> tau;
auto f = binFreq(n);
std::vector<double> e;
Complex nf3(0.0,0.0);
int j;
for ( j = 0 ; j < n ; j++ )
{
auto op = [](CArray::value_type v) {
return std::conj(v);
};
auto nf1 = ((x * x.apply(op)) * (y * y.apply(op)));
nf3 += nf1[j];
}
auto nf2 =std::sqrt(nf3);
auto nf =nf2/(double)n;
cout << "nf" << nf <<endl;
}
my main.cpp:
#include <complex>
#include <iostream>
#include <valarray>
#include <malloc.h>
#include <string>
#include <stdlib.h>
#include <fstream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <cmath>
#include <vector>
#include "fft.cpp"
#include "cs_delay.cpp"
#include "iterDelayEst.cpp"
using namespace std;
char filename[] = "myfile.txt";
char filename2[] = "myfile2.txt";
typedef std::complex<double> Complex;
typedef std::valarray <Complex> CArray;
/***********************************************************************************************
* function declarations
************************** ******************************************************** ***********/
void fft(CArray& x);
void ifft(CArray& x);
std::vector<double> binFreq(int n);
void cs_delay(CArray& x, int rate_hz, int delay_s, int n);
double iterDelayEst(int n, CArray& x, CArray& x2);
int main()
{
int dTest_samples;
int cTest;
int n=299;
int i;
int j;
double x [n];
double y [n];
int rate_hz=1;
int delay_s=30;
/*****************************getting x*******************************/
string line;
double Result;
ifstream myfile (filename);
if (myfile.is_open())
{
for ( i = 0 ; (i < n) && (myfile >> x[i]) ; ++i)
cout << line << 'n';
stringstream convert(line);
if ( !(convert >> Result) )
Result = 0;
x[i]=Result;
}
else cout << "Unable to open file";
/*****************************getting y******************************/
string line2;
double Result2;
ifstream myfile2 (filename2);
if (myfile2.is_open())
{
for ( i = 0 ; (i < n) && (myfile2 >> y[i]) ; ++i)
cout << line2 << 'n';
stringstream convert(line2);
if ( !(convert >> Result2) )
Result2 = 0;
y[i]=Result2;
}
else cout << "Unable to open file2";
/***********************************************************************/
/*********************for x******************/
Complex test[n];
for ( i = 0 ; i < n ; ++i )
test[i] = x[i];
CArray data(test,n);
/*********************for y******************/
Complex test2[n];
for ( j = 0 ; j < n ; ++j )
test2[j] = y[j];
CArray data2(test2,n);
// forward fft
fft(data);
std::cout << "fft" << std::endl;
for (int i = 0; i <n; ++i)
{
cout << data[i] << endl;
}
// inverse fft
ifft(data);
std::cout << std::endl << "ifft" << std::endl;
for (int i = 0; i <n; ++i)
{
std::cout << data[i] << std::endl;
}
cs_delay(data, 1, 6, n);
for (int i = 0; i <n; ++i)
{
std::cout << data[i] << std::endl;
}
iterDelayEst(n, data, data2 );
return 0;
}
the fft function:
//fft.cpp
using namespace std;
const double PI = 3.141592653589793238460;
/************************************************************************************************
* Cooley–Tukey FFT (in-place, divide-and-conquer) *
* Higher memory requirements and redundancy although more intuitive *
************************************************************************************************/
typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;
// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x)
{
const size_t N = x.size();
if (N <= 1) return;
// divide
CArray even = x[std::slice(0, N/2, 2)];
CArray odd = x[std::slice(1, N/2, 2)];
// conquer
fft(even);
fft(odd);
// combine
for (size_t k = 0; k < N/2; ++k)
{
Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
x[k ] = even[k] + t;
x[k+N/2] = even[k] - t;
}
}
result of my fft in matlab:
>> fft([0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0]')
9.0000 + 0.0000i
-3.3099 - 1.8957i
0.0000 - 5.0273i
-3.7722 + 5.1864i
1.0000 + 0.0000i
2.4654 + 1.0512i
-0.0000 - 1.4966i
-0.9547 - 0.4595i
1.0000 + 0.0000i
-1.3771 + 0.0371i
0.0000 - 0.6682i
-0.4381 + 1.8523i
1.0000 + 0.0000i
0.5733 - 0.8409i
0.0000 - 0.1989i
-1.1867 - 0.2275i
1.0000 + 0.0000i
-1.1867 + 0.2275i
0.0000 + 0.1989i
0.5733 + 0.8409i
1.0000 + 0.0000i
-0.4381 - 1.8523i
0.0000 + 0.6682i
-1.3771 - 0.0371i
1.0000 + 0.0000i
-0.9547 + 0.4595i
-0.0000 + 1.4966i
2.4654 - 1.0512i
1.0000 + 0.0000i
-3.7722 - 5.1864i
0.0000 + 5.0273i
-3.3099 + 1.8957i
result of my fft in c++:
(9,0)
(-3.3099,-1.89568)
(1.72253e-16,-5.02734)
(-3.77219,5.1864)
(1,0)
(2.46544,1.05123)
(2.94822e-16,-1.49661)
(-0.954746,-0.459467)
(1,0)
(-1.37707,0.0371386)
(-6.28295e-17,-0.668179)
(-0.438105,1.85232)
(1,0)
(0.573279,-0.840935)
(2.78992e-16,-0.198912)
(-1.18671,-0.227505)
(1,0)
(-1.18671,0.227505)
(1.72253e-16,0.198912)
(0.573279,0.840935)
(1,0)
(-0.438105,-1.85232)
(-3.82452e-17,0.668179)
(-1.37707,-0.0371386)
(1,0)
(-0.954746,0.459467)
(-3.67545e-17,1.49661)
(2.46544,-1.05123)
(1,0)
(-3.77219,-5.1864)
(-7.8049e-16,5.02734)
(-3.3099,1.89568)
Let assume the result of nf is correct, I added new lines new added lines :
for ( j = 0 ; j < n ; j++ ){
auto xcorr2=(fd2Tau * x.apply(std::conj));
auto xcorr3=(std::abs(xcorr2[i]))/nf;
cout << "xcorr3" << xcorr3[i] << endl;
}
unfortunately now I am receiving following errors messages:
In file included from myfft.cpp:16:0:
iterDelayEst.cpp: In function ‘double iterDelayEst(int, CArray&, CArray&)’:
iterDelayEst.cpp:61:30: error: no match for ‘operator[]’ (operand types are ‘std::complex<double>’ and ‘int’)
cout << "xcorr3" << xcorr3[i] << endl;
^
Can someone help me to sort out what I am doing wrong ?
Aucun commentaire:
Enregistrer un commentaire