mock-gr/algorithms/fun8_3_jamming_H.cpp

//
// File: jamming_H.cpp
//
// MATLAB Coder version            : 5.2
// C/C++ source code generated on  : 10-Mar-2023 15:16:29
//

// Include Files
#include "fun8_3_jamming_H.h"
#include "fun8_3_ifft.h"
#include "fun8_3_jamming_H_data.h"
#include "fun8_3_randn.h"
#include "rt_nonfinite.h"
#include "coder_array.h"
#include <cmath>
#include <iostream>
#include <iomanip>
#include "fun8_3_eml_rand_mt19937ar_stateful.h"



void jamming_H_initialize()
{
    eml_rand_mt19937ar_stateful_init_8_3();
    isInitialized_jamming_H = true;
}

// Function Declarations
static double rt_hypotd_snf(double u0, double u1);

// Function Definitions
//
// Arguments    : double u0
//                double u1
// Return Type  : double
//
static double rt_hypotd_snf(double u0, double u1)
{
  double a;
  double y;
  a = std::abs(u0);
  y = std::abs(u1);
  if (a < y) {
    a /= y;
    y *= std::sqrt(a * a + 1.0);
  } else if (a > y) {
    y /= a;
    y = a * std::sqrt(y * y + 1.0);
  } else if (!std::isnan(y)) {
    y = a * 1.4142135623730951;
  }
  return y;
}

//
// 生成射频噪声干扰信号
//  参数:
//  N - 信号长度
//  f_s - 系统采样率
//  A - 信号幅度
//  B_n - 射频信号要产生的带宽
//  f_0 - 调制信号频率
//  返回值:
//  s - 生成的射频噪声干扰信号
//
// Arguments    : double N
//                double f_s
//                double A
//                double B_n
//                double f_0
//                coder::array<creal_T, 2U> &s
// Return Type  : void
//
void jamming_H(double N, double f_s, double A, double B_n, double f_0,
               coder::array<creal_T, 2U> &s)
{
  coder::array<creal_T, 2U> S;
  coder::array<creal_T, 2U> s_n;
  coder::array<double, 2U> b_r;
  coder::array<double, 2U> c_b;
  coder::array<double, 2U> n;
  coder::array<double, 2U> r;
  coder::array<double, 1U> absdiff;
  double b_b[2];
  double N_n;
  double b;
  double bsum_im;
  double bsum_re;
  double xbar_im;
  int b_n;
  int bsum_re_tmp;
  int firstBlockLength;
  int k;
  int xblockoffset;
  N = 1024;
  f_s = 10e6;
  A = 1;
  B_n = 2e6;
  f_0 = 2e6;
  if (!isInitialized_jamming_H) {
    jamming_H_initialize();
  }
  if (std::isnan(N)) {
    n.set_size(1, 1);
    n[0] = rtNaN;
  } else if (N < 1.0) {
    n.set_size(1, 0);
  } else if (std::isinf(N) && (1.0 == N)) {
    n.set_size(1, 1);
    n[0] = rtNaN;
  } else {
    firstBlockLength = static_cast<int>(std::floor(N - 1.0));
    n.set_size(1, firstBlockLength + 1);
    for (bsum_re_tmp = 0; bsum_re_tmp <= firstBlockLength; bsum_re_tmp++) {
      n[bsum_re_tmp] = static_cast<double>(bsum_re_tmp) + 1.0;
    }
  }
  b = 1.0 / f_s;
  //  调制噪声参数
  // delta_F = B_n / 2; % 射频信号要产生的带宽B_n时，噪声要产生的带宽
  N_n = std::round(N * (B_n / (2.0 * f_s)));
  //  频谱上采样点数
  b_b[0] = 1.0;
  b_b[1] = static_cast<int>(N);
  coder::randn(b_b, r);
  firstBlockLength = static_cast<int>(N);
  b_r.set_size(1, static_cast<int>(N));
  for (bsum_re_tmp = 0; bsum_re_tmp < firstBlockLength; bsum_re_tmp++) {
    b_r[bsum_re_tmp] = r[bsum_re_tmp];
  }
  b_b[0] = 1.0;
  b_b[1] = static_cast<int>(N);
  coder::randn(b_b, r);
  firstBlockLength = static_cast<int>(N);
  c_b.set_size(1, static_cast<int>(N));
  for (bsum_re_tmp = 0; bsum_re_tmp < firstBlockLength; bsum_re_tmp++) {
    c_b[bsum_re_tmp] = r[bsum_re_tmp];
  }
  S.set_size(1, b_r.size(1));
  firstBlockLength = b_r.size(1);
  for (bsum_re_tmp = 0; bsum_re_tmp < firstBlockLength; bsum_re_tmp++) {
    S[bsum_re_tmp].re = b_r[bsum_re_tmp] + 0.0 * c_b[bsum_re_tmp];
    S[bsum_re_tmp].im = c_b[bsum_re_tmp];
  }
  //  高斯白噪声频谱,中值0，标准差1，维度(1，N)
  bsum_re_tmp = static_cast<int>(((N - N_n) - 1.0) + (1.0 - N_n));
  for (xblockoffset = 0; xblockoffset < bsum_re_tmp; xblockoffset++) {
    firstBlockLength =
        static_cast<int>(N_n + static_cast<double>(xblockoffset)) - 1;
    S[firstBlockLength].re = 0.0;
    S[firstBlockLength].im = 0.0;
  }
    coder::ifft_8_3(S, s_n);
  //  对做逆傅里叶变换就可以得到时域的有一定带宽的噪声信号
  b_n = s_n.size(1);
  if (s_n.size(1) == 0) {
    xbar_im = rtNaN;
  } else if (s_n.size(1) == 1) {
    if ((!std::isinf(s_n[0].re)) && (!std::isinf(s_n[0].im)) &&
        ((!std::isnan(s_n[0].re)) && (!std::isnan(s_n[0].im)))) {
      xbar_im = 0.0;
    } else {
      xbar_im = rtNaN;
    }
  } else {
    int lastBlockLength;
    int nblocks;
    if (s_n.size(1) <= 1024) {
      firstBlockLength = s_n.size(1);
      lastBlockLength = 0;
      nblocks = 1;
    } else {
      firstBlockLength = 1024;
      nblocks = s_n.size(1) / 1024;
      lastBlockLength = s_n.size(1) - (nblocks << 10);
      if (lastBlockLength > 0) {
        nblocks++;
      } else {
        lastBlockLength = 1024;
      }
    }
    N_n = s_n[0].re;
    xbar_im = s_n[0].im;
    for (k = 2; k <= firstBlockLength; k++) {
      N_n += s_n[k - 1].re;
      xbar_im += s_n[k - 1].im;
    }
    for (int ib{2}; ib <= nblocks; ib++) {
      xblockoffset = (ib - 1) << 10;
      bsum_re = s_n[xblockoffset].re;
      bsum_im = s_n[xblockoffset].im;
      if (ib == nblocks) {
        firstBlockLength = lastBlockLength;
      } else {
        firstBlockLength = 1024;
      }
      for (k = 2; k <= firstBlockLength; k++) {
        bsum_re_tmp = (xblockoffset + k) - 1;
        bsum_re += s_n[bsum_re_tmp].re;
        bsum_im += s_n[bsum_re_tmp].im;
      }
      N_n += bsum_re;
      xbar_im += bsum_im;
    }
    if (xbar_im == 0.0) {
      bsum_im = N_n / static_cast<double>(s_n.size(1));
      N_n = 0.0;
    } else if (N_n == 0.0) {
      bsum_im = 0.0;
      N_n = xbar_im / static_cast<double>(s_n.size(1));
    } else {
      bsum_im = N_n / static_cast<double>(s_n.size(1));
      N_n = xbar_im / static_cast<double>(s_n.size(1));
    }
    absdiff.set_size(s_n.size(1));
    for (k = 0; k < b_n; k++) {
      absdiff[k] = rt_hypotd_snf(s_n[k].re - bsum_im, s_n[k].im - N_n);
    }
    xbar_im = 0.0;
    N_n = 3.3121686421112381E-170;
    firstBlockLength = s_n.size(1);
    for (k = 0; k < firstBlockLength; k++) {
      if (absdiff[k] > N_n) {
        bsum_re = N_n / absdiff[k];
        xbar_im = xbar_im * bsum_re * bsum_re + 1.0;
        N_n = absdiff[k];
      } else {
        bsum_re = absdiff[k] / N_n;
        xbar_im += bsum_re * bsum_re;
      }
    }
    xbar_im = N_n * std::sqrt(xbar_im);
    xbar_im /= std::sqrt(static_cast<double>(s_n.size(1)) - 1.0);
  }
  s_n.set_size(1, s_n.size(1));
  firstBlockLength = s_n.size(1) - 1;
  for (bsum_re_tmp = 0; bsum_re_tmp <= firstBlockLength; bsum_re_tmp++) {
    N_n = s_n[bsum_re_tmp].re;
    bsum_re = s_n[bsum_re_tmp].im;
    if (bsum_re == 0.0) {
      bsum_im = N_n / xbar_im;
      N_n = 0.0;
    } else if (N_n == 0.0) {
      bsum_im = 0.0;
      N_n = bsum_re / xbar_im;
    } else {
      bsum_im = N_n / xbar_im;
      N_n = bsum_re / xbar_im;
    }
    s_n[bsum_re_tmp].re = bsum_im;
    s_n[bsum_re_tmp].im = N_n;
  }
  //  归一化,方差为1
  //  产生需要调制的射频噪声干扰s_0
  bsum_re = f_0 * 0.0;
  bsum_im = f_0 * 6.2831853071795862;
  s.set_size(1, n.size(1));
  firstBlockLength = n.size(1);
  for (bsum_re_tmp = 0; bsum_re_tmp < firstBlockLength; bsum_re_tmp++) {
    s[bsum_re_tmp].re = b * (n[bsum_re_tmp] * bsum_re);
    s[bsum_re_tmp].im = b * (n[bsum_re_tmp] * bsum_im);
  }
  firstBlockLength = s.size(1);
  for (k = 0; k < firstBlockLength; k++) {
    if (s[k].im == 0.0) {
      s[k].re = std::exp(s[k].re);
      s[k].im = 0.0;
    } else if (std::isinf(s[k].im) && std::isinf(s[k].re) && (s[k].re < 0.0)) {
      s[k].re = 0.0;
      s[k].im = 0.0;
    } else {
      N_n = std::exp(s[k].re / 2.0);
      s[k].re = N_n * (N_n * std::cos(s[k].im));
      s[k].im = N_n * (N_n * std::sin(s[k].im));
    }
  }
  s.set_size(1, s.size(1));
  firstBlockLength = s.size(1) - 1;
  for (bsum_re_tmp = 0; bsum_re_tmp <= firstBlockLength; bsum_re_tmp++) {
    N_n = A * s[bsum_re_tmp].re;
    bsum_re = A * s[bsum_re_tmp].im;
    s[bsum_re_tmp].re =
        N_n * s_n[bsum_re_tmp].re - bsum_re * s_n[bsum_re_tmp].im;
    s[bsum_re_tmp].im =
        N_n * s_n[bsum_re_tmp].im + bsum_re * s_n[bsum_re_tmp].re;
//	printf("%f  %f  ", s[bsum_re_tmp].re, s[bsum_re_tmp].im);
//	printf("%f\n", s[bsum_re_tmp]);
  }
}



//int main()
//{
//    double N=1024;
//    double f_s=10e6;
//    double A=1;
//    double B_n=2e6;
//    double f_0=2e6;
//    coder::array<creal_T, 2U> s;
//    jamming_H(N, f_s, A, B_n, f_0, s);
//    for (int i = 0; i < N; i++) {
//        std::cout << std::showpos << std::fixed << std::setprecision(6) << s[i].re << std::noshowpos << std::showpos << std::fixed << std::setprecision(6) << s[i].im << "i, ";
//    }
//    std::cout << "\n";
//    return 0;
//}