質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.78%

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

解決済

回答 1

投稿

  • 評価
  • クリップ 0
  • VIEW 858

aqufiz

score 65

前提・実現したいこと

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;
}
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • aqufiz

    2019/05/08 01:02 編集

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

    キャンセル

  • aqufiz

    2019/05/08 01:28

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

    キャンセル

  • episteme

    2019/05/08 04:31

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

    キャンセル

回答 1

checkベストアンサー

+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をチェックしてみるのが良いと思います。

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2019/05/08 12:14

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.78%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る

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