Search code examples
c++complex-numbersfrequency-analysis

C++ - Complex Value Calculation Error in Cross Spectral Density


I am a very beginner in c++ and I want to do some spectral calculations, in this case calculating the 'Cross Spectral Density' of two signal (vecFirst, vecSecond), which are already processed with a FastFourierTransformation. Resulting in freqvec and freqvec2, containing complex values for each frequency.

For this calculating it is essential to keep every value as a complex value. E.g.: CoSpectrum, which is calculated in line 6, should has a complex value as a result.

RowVectorXcd freqvec;
RowVectorXcd freqvec2;
fft.fwd(freqvec, vecFirst);
fft.fwd(freqvec2, vecSecond);

// # Create conjugate complex
freqvec.conjugate();
freqvec2.conjugate();
RowVectorXcd Rxy(freqvec.cols());
for (int i = 0; i < freqvec.cols(); i++) {
        std::complex<double>CoSpectrum( freqvec(i).real() * freqvec2(i).real() + freqvec(i).imag() * freqvec2(i).imag()) ;
        std::complex<double>QuadSpectrum( freqvec(i).real() * freqvec2(i).imag() - freqvec(i).real() * freqvec2(i).imag() ) ;
        std::complex<double>CoSpectrum_sqr = CoSpectrum * CoSpectrum ;
        std::complex<double>QuadSpectrum_sqr = QuadSpectrum * QuadSpectrum ;
        Rxy(i) = sqrt(std::complex<double>(CoSpectrum_sqr + QuadSpectrum_sqr)) ;
    }
}

Unfortunately I only get complex values with zero in the imaginary part.

Can anyone tell me why? I am guessing the expression freqvec(i).real() only returns a double value, but how can I get the real part but keep it a complexvalue. Or, accordingly, just multiply the imaginary part of a complex number with the real part of another and keep it the result a complex double.

Thanks for any help in advance.


Solution

  • At first: Thank you very much. I had thought of something like that, but wasn't sure. Thank got there a nice people like you who care about beginners.

    So I changed the snippet as following:

        // ### Attempting to compute the Frequency Power for Frequency Bins..
        RowVectorXcd freqvec;
        RowVectorXcd freqvec2;
        fft.fwd(freqvec, vecFirst);
        fft.fwd(freqvec2, vecSecond);
        std::cout<<"freqvec:"<<freqvec.cols()<<std::endl;
    
        // ### Attempting to compute the PowerSpectralDensitiy(PSD) and CrossSpectralDensity(CSD). The cross-spectral density is the Fourier transform of the cross-correlation function.
        RowVectorXcd Rxy(n_Epochs, freqvec.cols());
        RowVectorXcd Rxx(n_Epochs, freqvec.cols());
        RowVectorXcd Ryy(n_Epochs, freqvec.cols());
        for (int i = 0; i < n_Epochs; i++) {
                std::complex<double>CoSpectrum( std::complex<double>(freqvec(i).real(),0) * std::complex<double>(freqvec2(i).real(),0) + std::complex<double>(0,freqvec(i).imag()) * std::complex<double>(0, freqvec2(i).imag()) ) ;
                std::complex<double>QuadSpectrum( std::complex<double>(freqvec(i).real(), 0) * std::complex<double>(0, freqvec2(i).imag()) - std::complex<double>(0,freqvec(i).imag()) * std::complex<double>(freqvec2(i).real(), 0) ) ;
                std::complex<double>CoSpectrum_sqr = CoSpectrum * CoSpectrum ;
                std::complex<double>QuadSpectrum_sqr = QuadSpectrum * QuadSpectrum ;
                Rxy(i) = sqrt(std::complex<double>(CoSpectrum_sqr + QuadSpectrum_sqr)) ;
                Rxx(i) = std::complex<double>(freqvec(i).real(), 0) * std::complex<double>(freqvec(i).real(), 0) + std::complex<double>(0, freqvec(i).imag()) * std::complex<double>(0, freqvec(i).imag()) ;
                Ryy(i) = std::complex<double>(freqvec2(i).real(), 0) * std::complex<double>(freqvec2(i).real(), 0) + std::complex<double>(0, freqvec2(i).imag()) * std::complex<double>(0, freqvec2(i).imag()) ;
            }
        }
    

    This solved my problem. Thanks again for the nice discussion.