1
$\begingroup$

I am able to use the DFT and IDFT formulae to get the DFT of a window and reconstruct it but only if the sidelobes is very low (or the main lobe is narrow) as shown in the image below: enter image description here

Why does this happen? This "phenomenon" gets worse with higher sidelobes (such as rectangular window). I am not using the FFT/IIFT algorithm. I am using the formulae given in the Wikipedia page

Edit: It also gets worse if we increase the number of window samples.

Here is the code I used (it is not exactly DFT but it is "Beamforming DFT"):

clc;clear;close all

dx = 0.2;
freq = 900e6;
c = 299792458;          % speed of light in m/s
lambda = c/freq;

theta = -90:0.05:90;
pi_dx_by_lambda = dx*pi/lambda;
u = pi_dx_by_lambda*sind(theta);

Nelements = 30;

AmpTest = chebwin1(Nelements,50); ones(Nelements,1);

W_f = DFT(AmpTest,Nelements,u);
W_f_dB = 10*log10(abs(W_f).^2);
norm_W_f_dB = W_f_dB - max(W_f_dB);

Amp_new = IDFT(W_f,Nelements,u);

figure;subplot(2,2,1);plot(AmpTest,'r','Displayname','Original window');hold on;plot(Amp_new,'b--','Displayname','Window from IDFT');
ylabel('Amplitude');xlabel('Samples');title('Window function with 50dB SLL');legend
subplot(2,2,2);plot(theta,norm_W_f_dB);ylabel('Magnitude (dB)');title('DFT of Window');xlabel('Bins');


AmpTest = chebwin1(Nelements,25); ones(Nelements,1);
W_f = DFT(AmpTest,Nelements,u);
W_f_dB = 10*log10(abs(W_f).^2);
norm_W_f_dB = W_f_dB - max(W_f_dB);

Amp_new = IDFT(W_f,Nelements,u);

subplot(2,2,3);plot(AmpTest,'r','Displayname','Original window');hold on;plot(Amp_new,'b--','Displayname','Window from IDFT');
ylabel('Amplitude');xlabel('Samples');title('Chebyshev window function with 25dB SLL');legend
subplot(2,2,4);plot(theta,norm_W_f_dB);ylabel('Magnitude (dB)');title('DFT of Window');xlabel('Bins');

function W_f = DFT(Amp,Nelements,u)
% Compute the spectral window
W_f = zeros(1,length(u));
for k = 1:length(u)
    for n = 0:Nelements - 1
        W_f(k) = W_f(k) + Amp(n+1)*exp(1i*2*u(k)*n);
    end
end
end

function Amp = IDFT(W_f,Nelements,u)
% Compute the weights from the spectral window
Amp = zeros(1,Nelements);
for n = 0:Nelements - 1
    for k = 1:length(u)
        Amp(n+1) = Amp(n+1) + W_f(k)*exp(-1i*2*u(k)*n);
    end
end
IWm = real(Amp);
Amp = IWm./max(IWm);
end
$\endgroup$
3
  • $\begingroup$ Should I share my code? $\endgroup$ Commented May 4 at 8:45
  • 1
    $\begingroup$ Yes please share your code. Whether you use the FFT/IFFT or do it "brute force" as the DFT and inverse DFT the results should be identical: the FFT is just an efficient algorithm to give you the (exact) DFT result, not an approximation. $\endgroup$ Commented May 4 at 11:51
  • $\begingroup$ @DanBoschen, I have updated my question to include the code. Thanks for the help $\endgroup$ Commented May 6 at 2:33

2 Answers 2

3
+50
$\begingroup$

I think you are maybe misunderstanding how the "beamforming DFT" is actually computed. The beamforming equation you reference is \begin{equation} X[\theta] = \sum_{n=0}^{N-1}x[n]e^{j\frac{2\pi}{\lambda}nd\sin\left(\theta\right)} \end{equation} How this is practically computed is by letting $k = N\frac{d\sin(\theta)}{\lambda}$. The equation then becomes \begin{equation} X[k] = \sum_{n=0}^{N-1}x[n]e^{j2\pi\frac{kn}{N}} \end{equation} which is the DFT. We can then say $X[\theta] = X\begin{bmatrix}\arcsin\left(\frac{k\lambda}{Nd}\right)\end{bmatrix}$. This utilizes the mapping from "digital frequency" to "beamforming angle" $\theta = \arcsin\left(\frac{k\lambda}{Nd}\right)$. So, all you are doing is computing a DFT and then mapping it to the appropriate angle.

For the case when $d < \frac{\lambda}{2}$, $\arcsin$ becomes undefined at the larger angles. In order to avoid this, we scale the inside of the $\arcsin$ argument by $\frac{2d}{\lambda} $ to take into account the array's geometry proportional to the wavelength, because the apparent sample rate changes as you steer off boresight. For more info on this, see my answer here. This gives $\theta = \arcsin\left(\frac{2k}{N}\right)$.

As for your code, dx as defined is aliased. Did you mean dx = 0.2*lambda? Also, DFTs (not computed through FFT) is best represented as a matrix operation. To do this, you have to zero pad your signal to be the same length as the number of frequency sample points. Otherwise, you will end up with an underdetermined system. When you then do the IDFT, you take the first Nelements samples to remove the zero-padding. I have updated the code to show this.

clc;clear;close all

freq = 900e6;
c = 299792458;          % speed of light in m/s
lambda = c/freq;
dx = 0.2*lambda;

theta = -90:0.05:90;
pi_dx_by_lambda = 2*dx*pi/lambda;
u = pi_dx_by_lambda*sind(theta);

k = (0:length(theta)-1) - (length(theta)-1)/2;

Nelements = 30;

AmpTest = chebwin(Nelements,50); ones(Nelements,1);

W_f = DFT(AmpTest,Nelements,u);
W_f_dB = 10*log10(abs(W_f).^2);
norm_W_f_dB = W_f_dB - max(W_f_dB);

Amp_new = abs(IDFT(W_f,Nelements,u)); % abs() to take care of rounding errors

figure;subplot(2,2,1);plot(AmpTest,'r','Displayname','Original window');hold on;plot(1:Nelements,Amp_new(1:Nelements),'b--','Displayname','Window from IDFT');
ylabel('Amplitude');xlabel('Samples');title('Window function with 50dB SLL');legend
subplot(2,2,2);plot(asind(k*2/length(k)),fftshift(norm_W_f_dB));ylabel('Magnitude (dB)');title('DFT of Window');xlabel('Angle (degrees)');


AmpTest = chebwin(Nelements,25); ones(Nelements,1);
W_f = DFT(AmpTest,Nelements,u);
W_f_dB = 10*log10(abs(W_f).^2);
norm_W_f_dB = W_f_dB - max(W_f_dB);

Amp_new = abs(IDFT(W_f,Nelements,u)); % abs() to take care of rounding errors

subplot(2,2,3);plot(AmpTest,'r','Displayname','Original window');hold on;plot(1:Nelements,Amp_new(1:Nelements),'b--','Displayname','Window from IDFT');
ylabel('Amplitude');xlabel('Samples');title('Chebyshev window function with 25dB SLL');legend
subplot(2,2,4);plot(asind(k*2/length(k)),fftshift(norm_W_f_dB));ylabel('Magnitude (dB)');title('DFT of Window');xlabel('Angle (degrees)');

function W_f = DFT(Amp,Nelements,u)
    sz = max(Nelements,length(u));
    A = dftmtx(sz)/sqrt(sz);
    temp = Amp;
    Amp = zeros(sz,1);
    Amp(1:length(temp)) = temp;
    W_f = (A*Amp).';
end

function Amp = IDFT(W_f,Nelements,u)
    sz = max(Nelements,length(u));
    A = conj(dftmtx(sz))/sqrt(sz);
    temp = W_f;
    W_f = zeros(sz,1);
    W_f(1:length(temp)) = temp;
    Amp = (A*W_f).';
end
$\endgroup$
5
  • $\begingroup$ dx (element spacing) is 0.2m. Yes, I realize it should be lesser ( ~0.16m) to avoid grating lobes. $\endgroup$ Commented May 8 at 4:37
  • $\begingroup$ what does zero padding in the context of beamforming mean? Does it mean turning off some antenna elements (that is, they have zero amplitude)? $\endgroup$ Commented May 8 at 4:37
  • $\begingroup$ But even though this fixes the problem, why does the SLL effect the ability to reconstruct the window? $\endgroup$ Commented May 8 at 4:39
  • 1
    $\begingroup$ @RajaKrishnappa zero-padding would be the equivalent of imagining you had a lot more antenna elements but they are all set to 0. What it does is interpolate the beam pattern. If you have $M$ elements and $N$ beam pattern sample points, you violate the conditions of the DFT. You have to have $M=N$. The way to do this is zero-padding which interpolates the values appropriately. When you IDFT this sequence, you will get your original samples followed by a bunch of zeros. So, you just have to extract the correct sequence. $\endgroup$
    – Baddioes
    Commented May 8 at 5:09
  • 1
    $\begingroup$ @RajaKrishnappa as to why in the non-working case the SLL affected the response so much, it's hard to say. If you look at the 50 dB SLL case, the edges aren't quite right. My guess is that the lower sidelobes reduced the bias that you had from an incorrect transform. The higher sidelobe case increased the bias, making it more varied. But I don't know if that's right, just a guess. $\endgroup$
    – Baddioes
    Commented May 8 at 5:11
0
$\begingroup$

You are using a formulation that I am not familiar with, which I can see has several problems, only one of which you've recognised.

The problem you've noticed occurs because the standard discrete Fourier transform formulas that you see published around the place that you have used only work because the frequencies used result in an orthogonal basis set. In contrast, your set of frequencies fail to achieve this, and so your IDFT function does not invert your DFT function. To obtain a solution for a non-orthogonal basis set, you instead need to use a more general linear algebra approach, which can be achieved, at minimum, by replacing your IDFT function with the following lines of Matlab code:

A = exp(1i*2*(0:Nelements - 1)'*u);
Amp_new = real(W_f/A);  % the imaginary component is only non-zero due to numerical errors

However, several other problems I noticed are:

  1. You are wanting to fit 3601 frequencies to 30 data points, which is an underdetermined system of equations, hence has an infinite number of solutions. You should fit the same number of frequencies as data points (you can fit less frequencies, but then it becomes a least squares problem and you won't recover your signal when you transform then inverse transform).
  2. You are effectively doing the inverse transform from the time domain, and the transform from the frequency domain, which I maintained consistency with in my above Matlab code, but doesn't make any sense as at no point are you fitting frequency components to your data. This can't be fixed until problem (1) is fixed, that is, until it is no longer underdetermined.
  3. Your frequencies do not complete integer number of cycles in the window width + 1 period, which means, even though the mean of results may be correct, there will be a standard error.
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.