C++のfftを用いたofdmのプログラムで2の累乗でないものを0つめして2の累乗にしてから計算するとセグメントエラーが出る。

Question

前提・実現したいこと

C++のfftを用いたofdmのプログラムで2の累乗でないものを0つめして2の累乗にしてから計算するとセグメントエラーが発生しました。
入力データが元から2の累乗だとエラーは出ません。
0つめはqpskで変調した後、IFFTの直前で行なっています。

発生している問題・エラーメッセージ

セグメントエラー

試したこと

0つめが2累乗の個数までできていないのかと思い、確認しましたができていました。

該当のソースコード

#include<iostream>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<fstream>
#include<complex>
#include<algorithm>
#include<vector>

using namespace std;

#define PI acos(double(-1.0))
#define L 1024	/*送信データbit数*/
#define SNR_STEP 40/*SN比をいくつまで調べるか*/


typedef struct{
	vector<int> dat;
	vector<complex<double> > sig;
}TX_SIGNALS;

typedef struct{
	vector<int> dat;
	vector<complex<double> > sig;
	vector<complex<double> > noise;
}RX_SIGNALS;



void _randData(TX_SIGNALS *d){
	int i;
	double x;
	for(i=0;i<L;i++){
		x = drand48();			/*ランダムデータ発生*/
		if(x >=0.5)
			d->dat[i] =1.0;
		else
			d->dat[i] =0.0;
	}
}

void _qpskMod(TX_SIGNALS *s){
    int i,co;
    s->sig.resize(L/2);
    double c=1.0/sqrt(2.0);
    if ((s->sig.size()*2) != (s->dat.size())){
        cout<<"Error: _qpskMod, 長さが一致しません"<<endl;            /*QPSK変調器(送信)*/
        exit(1);
    }
    else {
        for(i=0;i<s->sig.size();i++){
            co = 2*i;
            complex<double>z(c * ((s->dat[co])-0.5)*2.0,c * ((s->dat[co+1])-0.5)*2.0);
            s->sig[i] =z;
        }
    }
}

int Ifour(TX_SIGNALS *tx,int isign){
    unsigned long n,mmax,m,j,istep,i,res,exs;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
    
    res=1;
    exs=L;
    if ((exs & (exs - 1)) != 0) {
        while (exs > 0) {
            res <<= 1; exs >>= 1;
        }
        tx->sig.resize(res/2);
        for (i = 0; i < (res - L)/2; i++) {
            const complex<double> add(0.0, 0.0);
            tx->sig.push_back(add);
        }
    }
    else{
        res=L;
    }
    n = res/2;
    
    j=0;
    for(i=0;i<n;i++){
        if(j>i){
            swap(tx->sig[j],tx->sig[i]);
        }
        m=n >>1;
        while (m >= 1 && j >=m){
            j -=m;
            m >>=1;
        }
        j+= m;
    }
    mmax=1;
    while(n >mmax){
        istep=mmax <<1;
        theta=isign*(2*PI/istep);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        for(m=0;m<mmax;m++){
            for(i=m;i<=n;i+=istep){
                j=i+mmax;
                tempr=wr*(tx->sig[j].real())-wi*(tx->sig[j].imag());
                tempi=wr*(tx->sig[j].imag())+wi*(tx->sig[j].real());
                complex<double> a((tx->sig[i].real())-tempr,(tx->sig[i].imag())-tempi);
                tx->sig[j]=a;
                complex<double> b((tx->sig[i].real())+tempr,(tx->sig[i].imag())+tempi);
                tx->sig[i]=b;
                
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }
    for(i=0;i<tx->sig.size();i++){
        tx->sig[i]/=(double)(tx->sig.size());
    }
    
    /*fstream file3;
    file3.open("ofIFFT.txt",ios::out);
    for(i=0;i<tx->sig.size();i++){
        file3<<(tx->sig[i]).real()<<"       "<<(tx->sig[i]).imag()<<endl;
        
    }
    file3.close();*/
    
    return(res);
    
}


double _SNRdBNP(double c_dB){
    return pow(10,(-1)*c_dB/10);                /*雑音電力計算*/
}


void _awgn(RX_SIGNALS *n,TX_SIGNALS *m, double Pn,double ave,double sigma){
    double B1,B2,B;
    int i;
    n->noise.resize(m->sig.size());
    for(i=0;i<(n->noise.size());i++){
        
        B1 =ave+ sigma*sqrt(-2.0*log(drand48()))*cos(2.0*PI*drand48()); /*Box Muller法*/
        B2 =ave+ sigma*sqrt(-2.0*log(drand48()))*sin(2.0*PI*drand48());
        
        complex<double>B(B1* sqrt(Pn/2),B2* sqrt(Pn/2));    /*AWGN発生*/
    
        
        n->noise[i]=B;
    }
    /*fstream file4;
    file4.open("noise.txt",ios::out);
    for(i=0;i<n->noise.size();i++){
        file4<<n->noise[i]<<endl;
    }
    
    file4.close();*/
}

void _vectorSum(RX_SIGNALS *a,TX_SIGNALS *b){
    int i;
    a->sig.resize(b->sig.size());
    if(a->sig.size() != b->sig.size() || b->sig.size() != a->noise.size()){
        cout<<"Error: _vectorSum, 長さが一致しません"<<endl;    /*ベクトル加算（系列と系列の加算）*/
        exit(1);
    }
    else{
        for(i=0;i<(a->sig.size());i++){
            complex<double> H(1+pow(10,-0.25)*cos(2*PI*5*pow(10.0,-6.0)*i),-pow(10,-0.25)*sin(2*PI*5*pow(10.0,-6.0)*i));
            a->sig[i] = H*(b->sig[i]) + a->noise[i];
        }
    }
}

void four(RX_SIGNALS *rx,int isign,int nn){
    unsigned long n,mmax,m,j,istep,i,res,exs;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
    
    n=nn/2;
    
    j=0;
    for(i=0;i<n;i++){
        if(j>i){
            swap(rx->sig[j],rx->sig[i]);
        }
        m=n >>1;
        while (m >= 1 && j >=m){
            j -=m;
            m >>=1;
        }
        j+= m;
    }
    mmax=1;
    while(n >mmax){
        istep=mmax <<1;
        theta=isign*(2*PI/istep);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        for(m=0;m<mmax;m++){
            for(i=m;i<=n;i+=istep){
                j=i+mmax;
                tempr=wr*(rx->sig[j].real())-wi*(rx->sig[j].imag());
                tempi=wr*(rx->sig[j].imag())+wi*(rx->sig[j].real());
                complex<double> a((rx->sig[i].real())-tempr,(rx->sig[i].imag())-tempi);
                rx->sig[j]=a;
                complex<double> b((rx->sig[i].real())+tempr,(rx->sig[i].imag())+tempi);
                rx->sig[i]=b;
                
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }

}


void zero(RX_SIGNALS *d){
    int i;
    for(i=0;i<(d->sig.size());i++){
        complex<double> H(1.0+pow(10.0,-0.25)*cos(2.0*PI*5.0*0.000001*i),-pow(10.0,-0.25)*sin(2.0*PI*5.0*0.000001*i));
        d->sig[i]=d->sig[i]/H;
    }
}

void _qpskDem(RX_SIGNALS *rx){
    int i,co;
    for(i=0;i<rx->dat.size();i++)
        rx->dat[i] = 0.0;
    for(i=0;i<rx->sig.size();i++){
        co = 2*i;
        if((rx->sig[i]).real()>=0.0)                    /*QPSK復調器(受信)*/
            rx->dat[co]=1;
        if((rx->sig[i]).imag()>=0.0)
            rx->dat[co+1]=1;
    }
}


double _BER(TX_SIGNALS *tx,RX_SIGNALS *rx){
	int i,sum=0;
	double BER;
	if(tx->dat.size() != rx->dat.size()){
		cout<<"Error: _BER() 長さが一致しません"<<endl;	/*BER計算*/
		return 1;
	}
	else{
		for(i=0;i<tx->dat.size();i++)
			if((tx->dat[i]) - (rx->dat[i])!=0){
				sum++;
			}

		BER = (double)sum/(tx->dat.size());
	}
	return(BER);
}

void _BERprint(int snr_step, double *SNRdB, double *BER){
	int i;
	fstream file1;
	file1.open("OFDMBER.txt",ios::out);					/*BERをファイルに書き出す*/


		for(i=0;i<snr_step;i++){

			file1<<*(SNRdB+i)<<" "<<*(BER+i)<<endl;
		}
	file1.close();
}



int main(){
	TX_SIGNALS tx;
	RX_SIGNALS rx;
	double BER[SNR_STEP],SNRdB[SNR_STEP];
	int i,y,news;
	/*init_rtx(&tx,&rx);		TX_SIGNALS,RX_SIGNALSのメンバ初期化*/

	tx.dat.resize(L);
	tx.sig.resize(L/2);	/*txのメモリの確保*/
	rx.dat.resize(L);
	rx.sig.resize(L/2);
	rx.noise.resize(L/2);	/*rxのメモリの確保*/
	srand48(123);
	for(i=0;i<SNR_STEP;i++){
		SNRdB[i]=(double)i;
	}
	for(i=0;i<SNR_STEP;i++){
		_randData(&tx);		/*ランダムデータ発生*/
		_qpskMod(&tx);			/*QPSK変調器(送信)*/
		news=Ifour(&tx,-1);			/*Ifft*/
		_awgn(&rx,&tx,_SNRdBNP(SNRdB[i]),0.0,1.0);	/*AWGN発生*/
		_vectorSum(&rx,&tx);	/*ベクトル加算（系列と系列の加算）*/
		four(&rx,1,news);				/*fft*/       
		zero(&rx);					/*等化*/
		_qpskDem(&rx);			/*QPSK復調器(受信)*/
		BER[i]=_BER(&tx,&rx);			/*BER計算*/
	}
	_BERprint(SNR_STEP,SNRdB,BER);
	return 0;
}

Comments

[moredeep]

エラーが出た時にどの行でエラーが発生したか表示されていませんか？

[aqufiz]

すみません、セグメントエラーなので表示されてないです。

[aqufiz]

各関数が終わった時にコメントを表示させるようにした時、fftのところで止まっているので、fftのプログラムに問題があるのかもしれないです。

[moredeep]

すみませんがちょっとよくわからないんです。
・入力データとはどれのことを指していますか？
・0つめとは何のことですか？何を目的としてどのような処理を行おうとしていますか？
・提示のコードで動くことを確認しましたか？
→exs=L;←Lが定数(1024)なので、その直後のifは必ず偽です。
→その下のwhile内のifの条件(i<=n)、tx->sigのサイズとnはたぶん同じです。iがnと一致してしまうと、tx->sig[tx->sig.size()]と同じ意味になり、out of rangeです。iはm(初回は0)から、nは512、iStepはmmax(初回は1)+1なので、必ず1ループ目にエラーです。

他にもあれ？と思うようなところがあり、でも全てを見たわけではないのでもしかしたら正常かもしれないと、恐らく全体を把握しないと判断が難しいのですが、さすがにこのレベルのコード全部は見ていられないです。。

[moredeep]

drand48を使っているということはUnix系だと思いますので、以下のサイト辺りを参考に特定できませんかね？
http://rabbitfoot141.hatenablog.com/entry/2016/11/14/153101

[episteme]

> 入力データが元から2の累乗だとエラーは出ません。
0で埋めた直後、データ数が確かに2のベキになってるかは確認済なんですよね?

[aqufiz]

>>moredeepさん
・入力データは送信データビット数のLを示しています．
・fftは2の累乗でないと計算できないため，足りないデータ数の分のところに0を代入して数合わせをしています
・Lは送信ビット数でいろいろな値に変えたりしています．
・i=nのときout of rangeになってしまっていました，ありがとうございます．
説明不足ですみません．お忙しい中，助言をしてくださってありがとうございます．

[aqufiz]

>>episteme
確認済みです．

[episteme]

↑だとするとデータ数がハナっから2^Nなときと0詰めしたときとでは何かが違うってことやね...

[aqufiz]

そうですね、データに0が入っていること以外同じはずなのですが、なぜかセグメントエラーが出てしまいます。

[aqufiz]

i=nのときout of rangeになってしまっていたところをIFFTだけでなく，FFTでも修正したところコンパイルが通り，結果も直っていましたありがとうございます．
moredeepさんを回答してくれた方にしたいので，回答の方にコメントしていただけると助かります．

[episteme]

「ハナっから2^Nならちゃんと動く」を説明できるんですかその↑修正?

[aqufiz]

ご指摘の通りだと思います。一応動いていますが求めてる動作ではないです。しかし、今回の質問ではセグメントエラーの解決を目指していたので解決としました。

[aqufiz]

2^Nの時も求めた結果ではなかったのですが、2^Nでないデータを入れた時に結果すら見れなくなってしまったので、まずそこを修正してからより深い、白色雑音の大きさなどの中身の修正をしていきたいと思い今回の質問をさせていただきました。

[aqufiz]

こんな情報も少なく曖昧な質問に答えていただきありがとうございました。

[episteme]

fftを実装することが目的でないのなら出来合いのFFTWあたりを使っちゃうのもテかも。老婆心ながら。

Answer 1

本当に根本から解決できたかは置いておき、とりあえずエラーが出なくなったとのことで、お役に立てて良かったです。

他にも問題と思しき箇所、例えば

    if ((exs & (exs - 1)) != 0) {
        while (exs > 0) {
            res <<= 1; exs >>= 1;
        }
        tx->sig.resize(res/2); // ←これ
        for (i = 0; i < (res - L)/2; i++) {
            const complex<double> add(0.0, 0.0);
            tx->sig.push_back(add); // ←これ
        }
    }

入力Lが6のとき、whileの後resは8であると思います。
tx->sigのサイズはres / 2で4ですよね。
その後に(res - L) / 2回((8-6)/2で1回)push_backすると、最終的にはtx->sigのサイズは5です。これは意図した動きではないのではないかと思います。
(どんな入力の時にこうなってほしいとかが曖昧なのでちょっと自信ないですが。)

兎に角、一度関数毎のIOをチェックしてみるのが良いと思います。

Comments

[aqufiz]

tx->sig.resize(res/2); は必要なかったので削除しました。
ここからは関数の動作を一つずつ確認していきたいと思います。
ありがとうございました。