C ++中数组的复数共轭转置

Question

我正在尝试编写与以下 Matlab 代码等效的 C++ 代码：

n = 32;
s1 = dlmread('ssa.txt');
s2 = dlmread('ssb.txt');
fd1 = fft(s1);
fd2 = fft(s2);
nf = sqrt((fd1 * fd1') * (fd2 * fd2')) / n

在哪里

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

和

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

不幸的是，我在 C++ 中得到的结果与 Matlab 有点不同

来自 Matlab 的结果：

nf = 12

C++ 中的结果：

nf(4.88621,0)

我的 C++ 代码：

//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;
}

我的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;
}

fft 函数：

//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;
    }
}

我在 matlab 中的 fft 结果：

>> 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

我在 C++ 中的 fft 结果：

(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)

假设 nf 的结果是正确的，我添加了新行新添加的行：

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;
}

不幸的是，我现在收到以下错误消息：

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;


                              ^

谁能帮我解决我做错了什么？

我更正了代码，我不再收到错误了。但是，我添加了新行，我想将此 Matlab 代码（ ix = find(xcorr == max(xcorr)); ）转换为 C++ 该代码查找对应于 xcorr 元素最大值的值的索引。为此，我编写了以下代码：

find.cpp:
#include <iostream>
#include <vector>
#include <cstdio>
#include <algorithm>
#include <iterator>
#include <cstdlib>
#include <fstream>

using namespace std;
int main() {
    std::vector<double> v;
    v.push_back(1);
    v.push_back(3);
    v.push_back(5);
    v.push_back(0);
    v.push_back(7);
    v.push_back(1);
    v.push_back(7);

    std::vector<double>::iterator result ;
    const double maxim = *max_element(v.begin(), v.end());
    for(result = std::find(v.begin(), v.end(), maxim); result != v.end(); result = std::find(result + 1, v.end(), maxim))
    {
           cout << "maximum found at index " << result -v.begin() +1 << endl;
    }
    return 0;
}

我需要将此代码合并到我的 iterDelayEst.cpp 函数中，我不知道如何将 xcorr 的结果（数组）作为参数发送到 find 函数，以便它用作向量

Answer 1

正确代码：

double iterDelayEst(int n, CArray& x, CArray& y)
{

    //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 nfa=(x * x.apply(std::conj));
        nf4 +=nfa[j];

        auto nfb=(y * y.apply(std::conj));
        nf5 +=nfb[j];

    }

    auto nf1 = (nf4 * nf5);
    auto nf2 =std::sqrt(nf1);
    auto nf =nf2/(double)n;
    cout << "nf1" << nf1 <<endl;
    cout << "nf2" << nf2 <<endl;
    cout << "nf" << nf <<endl;
}

C ++中数组的复数共轭转置

The complex conjugate transpose of an array in c++

c++

matlab

c++14

transpose

complex-numbers