0

There are several thread regarding zero padding at the center and padding at the end.

I have a spectrum and want to interpolate it.

DFT of real sequence: $\left[X_0, X_1, X_2,\dots,-X_2, -X_1\right]$

Padding zeros: Now I can put zeros of required length from Nyquist till first -ve frequency bin.

My confusion is whether I should rearrange the spectrum like -ve,DC, +ve, Nyquist and padd at either side of the spectrum.

What is the correct way to do it?

int N                                                           = 9;

fftw_complex s_m_x_in                                          = (fftw_complex) fftw_malloc(sizeof(fftw_complex) * N);
fftw_complex s_m_x_out                                         = (fftw_complex) fftw_malloc(sizeof(fftw_complex) * N);


fftw_plan p1                                                    = fftw_plan_dft_1d(N, s_m_x_in, s_m_x_out, FFTW_FORWARD, FFTW_ESTIMATE);


s_m_x_in[0][0]                                                  = 2.0;
s_m_x_in[1][0]                                                  = 3.0;
s_m_x_in[2][0]                                                  = 6.0;
s_m_x_in[3][0]                                                  = 8.0;
s_m_x_in[4][0]                                                  = 1.0;
s_m_x_in[5][0]                                                  = 7.0;
s_m_x_in[6][0]                                                  = 3.0;
s_m_x_in[7][0]                                                  = 8.0;
s_m_x_in[8][0]                                                  = 8.0;


std::cout <<"==================original spectrum===================" <<std::endl;


fftw_execute(p1);


for (int i=0; i<N; i++)
{
    std::cout << s_m_x_out[i][0]<<"\t"<<s_m_x_out[i][1]<<std::endl;
}


std::cout << std::endl;


fftw_destroy_plan(p1);


int input_length_by_2                                           = (N % 2 == 0) ? N / 2 : (N-1) / 2;


int N_2                                                         = 2*N;


fftw_complex g_m_x_in                                          = (fftw_complex) fftw_malloc(sizeof(fftw_complex) * N_2);
fftw_complex g_m_x_out                                         = (fftw_complex) fftw_malloc(sizeof(fftw_complex) * N_2);


fftw_plan q1                                                    = fftw_plan_dft_1d(N_2, g_m_x_in, g_m_x_out, FFTW_BACKWARD, FFTW_ESTIMATE);


for (int i=0; i<=input_length_by_2; i++)
{
    g_m_x_in[i][0]                                              = s_m_x_out[i][0];
    g_m_x_in[i][1]                                              = s_m_x_out[i][1];


if (i &lt; input_length_by_2)
{
    g_m_x_in[(N_2 - 1) - i][0]                              = s_m_x_out[N - 1 - i][0];
    g_m_x_in[(N_2 - 1) - i][1]                              = s_m_x_out[N - 1 - i][1];
}                



}


std::cout <<"===========g_m_x_in (zeros at the center)=============" <<std::endl;


fftw_execute(q1);


for (int i=0; i<N_2; i++)
{
    std::cout << g_m_x_in[i][0]<<"\t"<<g_m_x_in[i][1]<<std::endl;
}


std::cout << std::endl;
std::cout <<"===========g_m_x_out (interpolated)===================" <<std::endl;


for (int i=0; i<N_2; i++)
{
    std::cout << g_m_x_out[i][0] / (N) <<"\t"<<g_m_x_out[i][1] / (N) <<std::endl;
}


std::cout <<"======================================================" <<std::endl;
std::cout << std::endl;


fftw_destroy_plan(q1);


fftw_complex c_m_x_in                                          = (fftw_complex) fftw_malloc(sizeof(fftw_complex) * N_2);
fftw_complex c_m_x_out                                         = (fftw_complex) fftw_malloc(sizeof(fftw_complex) * N_2);


fftw_plan q2                                                    = fftw_plan_dft_1d(N_2, c_m_x_in, c_m_x_out, FFTW_BACKWARD, FFTW_ESTIMATE);


int output_length_by_2                                          = (N_2 % 2 == 0) ? N_2 / 2 : (N_2 - 1);


for (int i=0; i<output_length_by_2; i++)
{
    if (i >= input_length_by_2)
    {
        c_m_x_in[output_length_by_2 - input_length_by_2 + i][0]     = s_m_x_out[i - input_length_by_2][0];
        c_m_x_in[output_length_by_2 - input_length_by_2 + i][1]     = s_m_x_out[i - input_length_by_2][1];

    }


if (i &lt; input_length_by_2)
{            
    c_m_x_in[output_length_by_2 - input_length_by_2 + i][0]     = s_m_x_out[i + input_length_by_2 + 1][0];
    c_m_x_in[output_length_by_2 - input_length_by_2 + i][1]     = s_m_x_out[i + input_length_by_2 + 1][1];

}                



}


std::cout <<"===========c_m_x_in (spectrum at the center)==========" <<std::endl;


fftw_execute(q2);


for (int i=0; i<N_2; i++)
{
    std::cout << c_m_x_in[i][0]<<"\t"<<c_m_x_in[i][1]<<std::endl;
}


std::cout << std::endl;
std::cout <<"===========c_m_x_out (interpolated)===================" <<std::endl;


for (int i=0; i<N_2; i++)
{
    std::cout << c_m_x_out[i][0] / (N) <<"\t"<<c_m_x_out[i][1] / (N) <<std::endl;
}


std::cout <<"======================================================" <<std::endl;
std::cout << std::endl;


fftw_destroy_plan(q2);


fftw_free(s_m_x_in);
fftw_free(s_m_x_out);


fftw_free(g_m_x_in);
fftw_free(g_m_x_out);


fftw_free(c_m_x_in);
fftw_free(c_m_x_out);
jomegaA
  • 659
  • 3
  • 16
  • If you want to interpolate the spectrum by zero-padding, the zero-padding is done in the time domain, not the frequency domain. Am I misunderstanding your intent? – Jdip Nov 27 '23 at 06:18
  • I think spectrum can be interpolated perhaps to get a nice graph in Time-Domain. – jomegaA Nov 27 '23 at 06:41
  • So you want to interpolate the time domain then? – Jdip Nov 27 '23 at 07:06
  • its the part of the algorithm where the samples of sphere from lower to higher level are interpolated. – jomegaA Nov 27 '23 at 07:22
  • updated the code snippet. Though it interpolates but when spectrum is symmetrical the interpolated values are negative compared to the zero-padding at the center. I do not know the reason... – jomegaA Nov 27 '23 at 10:03

0 Answers0