unwanted bandstop in kspaceLineRecon (problem in DCT computation)

Hi everybody,
if I am not mistaken, there is slight but significant mistake in the calculation of the DCT in the photoacoustic reconstruction script kspaceLineRecon.
The DCT is realized by temporally mirroring the measurement data to obtain an even signal and then performing an FFT. Therefore, a flipped version of the measurement data is concatenated with the actual measurement data while leaving out the first sample once:
p = [flipdim(p, 1); p(2:end, :)];
In this new dataset, the element corresponsing to t=0 is at (N+1)/2, if N is the (uneven) size of the new dataset. This operation is followed by an fftshift (before Fourier transform):
p = sf.*fftshift(fftn(fftshift(p,1)));
However, to my understanding, the fftshift here should be an ifftshift. The fftshift for an uneven number of samples (which is always the case here) returns the element with index (N+1)/2+1 as new first element, which corresponds to t=dt and not t=0. An ifftshift returns the element with index (N+1)/2 as new first element, which corresponds to t=0. For a Cosine transform to be computed as a Fourier transform, however, the signal has to be even, which implies that it has to be symmetric around the first element and thus the first element has to be that corresponding to t=0.
This does not sound like a big deal, but in fact, it leads to a situation, where the reconstructed signal after inverse Fourer transform is splitted up in imaginary part and real part and since the imaginary part is eliminated in the end, some of the actual content of the signal is lost, which can be seen when looking at the axial spectrum of reconstructed data of a broadband signal. There is a bandstop in the spectrum around fs/4, which means this frequency is not reconstructed at all. For a center frequency chosen as fs/4, which should not be too unusual, this implies a supproession of the most dominant frequnecyband.
To cut a long story short, I suggest to replace the inner fftshift by an ifftshift in:
p = sf.*fftshift(fftn(fftshift(p)));
such that:
p = sf.*fftshift(fftn(ifftshift(p)));
and consequently replace the outer ifftshift by an fftshift in:
p = real(ifftshift(ifftn(ifftshift(p))));
such that:
p = real(fftshift(ifftn(ifftshift(p))));
which should completely solve this problem.
I hope my explaination was not too confusing. Please let me know if I am misunderstanding anything here!
Thanks,
Hans-Martin

k-Wave

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields

unwanted bandstop in kspaceLineRecon (problem in DCT computation)

Reply

k-Wave

A MATLAB toolbox for the time-domain simulation of acoustic wave fields