Is it possible to 'interpolate' in FFT data?

We kmow OTF and PSF are continuous， can be interpolate inside it (<cutoff freq ) even when we express them wIth matrix use integer row&col coordinate .

now I have a image, and get its FFT2, and find some feature like max or min only with row and col,

but I know true max or min not at that integer row and col exactly， It is really at a fraction position near that integer row and col in matrix.

How to get it if I know fraction or decimal position ?

In frequency domain， can exten matrix size with 0, but can not interpolate although we can interpolate at time domain.

How to do ?

Maybe I can do so like such method : in time domain , matrix .* with exp( - j ke), {ke is a decimal number part }, to realize freq shift effect， after that, FFT, now at that row&col coordinate really is a value corresponding to [ row&col + decimal ke ] position, so I get it .

Such idea is right ? I am not expert at FFT, who can tell me , Is it Right ?

or how to find the feature not exactly at discrete integer row&col coordinate ?

I need some value at real flot number row&col coordinate ?

-----------------------------------------------------------------------------------------------------------------

I mean :

For simple example:

there are 512 data point with coordinate 1~512, transfer to freq domain with FFT, it is 512 mag and phase with coordinate -256~255 * 1/256 .

now if I specially add a cosine wave with freq 10.324 *1/512 to 512 data point (I used it as a special inner reference), then transfer to freq domain with FFT.

here , how to find mag and phase information of this added cosine wave in freq domain ?

freq 10.324 , is not a integer -256~255 coordinate, so how to find it ?

now sample code in 1D list below. ( I need work in 2D in real world) ： Is it possible? how?

% Time domain
t = linspace(1,512,512);
g = 2*sin(2*pi*t*16/512) + 0.5*sin(2*pi*t*64/512)+ 0.4*sin(2*pi*t*128/512);
% g = 2*sin(2*pi*t*16/512);
gadd = g+0.2*cos(2*pi*t*180.327/512);  
% gadd = g + sin(2*pi*t*10.327/512);  
% now add 180.324 or can be any another not N/512 

subplot(2,2,1);
plot(g);
title("Time domain 3 sin")
subplot(2,2,2);
plot(g);
title("Time domain, 3 sin + cos")

% freq domain
fco = linspace(-256,255,512);
F=fftshift(fft(g));
F1=fftshift(fft(gadd));
subplot(2,2,3);
stem (fco,real(F));
title("Freq domain, 3 sin")
subplot(2,2,4);
stem (fco,real(F1));
title("Freq domain, 3 sin + cos")

% how to find that 180.327  accuratly ?
% watch freq plot, this wave's energe distrbute on 180 and 181 two point,but how to get  180.327 's ?
% condition :  prior I know that freq wave inside them , how to get its information like mag or phase ?
%  Maybe should I use Fourier formula directly  by myself ？
%  like  result = sum(data.*exp(-i*2*pi*k*t));     
% k is I needed freq, only not any  n/datalength  limited by data count  --- integer numbers!

% if try ：
   result = sum(gadd.*exp(-1i*2*pi*t*180.327/512));     
% now get result is 50.2839072265070 - 0.424971233625617i
%  How  relate this complex number with wave's  Amp=0.2, and phase =0  ?
% try try
Amp=2*real(result)/512;
theta=angl(result);
% get 0.196 and -0.0085, compare with 0.2 and 0, 
% not too bad,  % I will try other.....  Maybe OK can be usded.  
% but  still  do not know how to do for 2D ?  my god !

[Rik]

It doesn't sound like it is really relevant that your data represents the result of an FFT. It sounds like you want to fit a surface to your data (at least locally) to find the location of a maximum on a sub-pixel level. Can you post the 2D array along with the code you tried so far?

[Adam]

I don't understand a lot of what you are asking, but yes, in theory you can interpolate in the frequency domain just as you can in the time domain. It comes with the same caveats, that you are effectively introducing data where there wasn't any before, using an algorithm of your choice (e.g. linear, pchip, spline, etc), but so long as the frequency space is relatively smooth (as has been the case when I have interpolated in it) it is no different to interpolating in any other space, such as the time domain.

I have interpolated in 2d in the joint time-frequency space in the past, to give 'more accurate' locations of the peaks in that space.

Whether it makes sense in your exact situation or not I couldn't really say, just that in a general case there's no more reason against it than for any other sufficiently smooth domain. Obviously if it isn't smooth then interpolating is always risky, because the true data could be well off what you interpolate.

[Mathieu NOE]

the result of fft2 can be viewed as an image so you could use this fex :

Fast 2D peak finder - File Exchange - MATLAB Central

with the sub pixel resolution option (weighted centroids)

[xd]

Fast 2D peak finder is in image.

I need is find in frequency domain of image, not in image. But in frequency domain, only in integer row&col coordinate. I need more accurate to find in frequency domain .

[Star Strider]

One approach may be to increase the number of frequencies the Fourier transform is evaluated at. The easiest way to do that is to zero-pad the original data. This will increase the frequency resolution of the Fourier transform in the relevant dimension (or dimensions), so you may be able to actually calculate the value you want, rather than needing to interpolate it. Both the fft and fft2 functions offer this option.

The frequency limits themselves will not change, and will of necessity be between zero and the Nyquist frequency in the relevant dimensions. You can restrict these frequency limits after you calculate the Fourier transform, however you cannot exceed them.

The greatest computational effiiciency will result if the size value for fft (or values for fft2), are integer powers of 2, (use the nextpow2 function to determine the closest value, although you may increase it), since that corresponds to the way the Fast Fourier Transform iis coded.

It will be necessary for you to experiment with this in order to get the result you want.

[Rik]

I would suggest reading Star Striders comment carefully, but if you still want to find the peak: a 2D Fourier transform of an image is also still a 2D array. If you are looking for a peak in a 2D array, you can use 2D array tools. Generally those tools will mention images, as that is the most common 2D array. However, the fact that your axes are frequency instead of spatial don't matter much for such functions.

[xd]

I am sorry. Now I have added sample code with 1D.

your opinion is based on textbook, but now my special situation is I know I added a standard reference wave to it, and I want to find it out more accuratly .

yes, there are limited by Nyquist rule, but in real world， theory is used as reference, and real problem is need to be resolved in a some accepted uncertainty.

And as I understand, this is limited by integer diacrete， It's not that it can't be done.

[Rik]

It is not clear to me who you're replying to and what your currently remaining issue is.

If you have a 2D problem, you should provide a 2D example. I stand by my opinion that you can just use the 2D peak finder. If you want more accuracy than that, you will have to look into what Star Strider told you and not dismiss it out of hand.

[xd]

my problem :

I added a standard wave to signal in time domain, want to find its feature in Freq domain.

But this wave frequency is not just at row&rol discrete position of DFT.

So I need to build a model to find it first from simple , but Discrete FT give it a dispersed expression.

So like that in code , I want to find the mag and phase of that cosine wave I added..

total experiment like this:

1) used some signal and my added cosine wave, to excitate a system, then record the reponse results

2) that added cosine wave used as a standard reference to quantify other things。

3） I know its freq not change , should find it out . but now, DFT calculation dispersed it .

[Mathieu NOE]

again , show us a working example if you have it

[Paul]

Essentially the zero padding approach shown here. Will become problematic if any other significant component of the signal has frequency too close to the frequency of the added cosine, in which case I'm not sure what could be done.

% Time domain
%t = linspace(1,512,512);
t = 0:511;
g = 2*sin(2*pi*t*16/512) + 0.5*sin(2*pi*t*64/512)+ 0.4*sin(2*pi*t*128/512);
% g = 2*sin(2*pi*t*16/512);
% add some noise
gadd = g + 0.2*cos(2*pi*t*180.327/512) + .2*randn(size(t));
[h,w] = freqz(gadd,1,2^14);
figure
plot(subplot(211),w*512/2/pi,real(h))
xline(180.327);
plot(subplot(212),w*512/2/pi,imag(h))
xline(180.327)

●

copyobj(get(gcf,'children'),figure);
set(get(gcf,'Children'),'XLim',[175,185])

[xd]

Thanks Paul ！ very good！

freqz(gadd,1,2^14);

is excellent ！

but maybe 2^14 need a lot calculation time.

I have no any knowledge about filter and signal spectrum analysis.

maybe i need to learn some.

Is it possible to calculate quickly ?

I will try to compare it with

sum(gadd.*exp(-1i*2*pi*t*180.327/512)); 

[Paul]

Originally, I thought the issue is that you don't really know the frequency of the added reference wave and that you're trying to find it through spectral analysis, hence the suggested computation of the DTFT of the signal using freqz, which could also be done with fft with zero-padding.

If you know that the reference frequency is 180.327/512 cycles/sample, then I suppose we are trying to estimate the reference wave's amplitude and phase. If that's the case, then your approach is reasonable, as long as the original signal is sufficiently separated in frequency from the reference and you correct for the amplitude and phase induced by the rectangular window at the reference frequency.

rng(100);
t = 0:511;
g = 2*sin(2*pi*t*16/512) + 0.5*sin(2*pi*t*64/512)+ 0.4*sin(2*pi*t*128/512);
% g = 2*sin(2*pi*t*16/512);
% add some noise
noise = 0.2*randn(size(t));
gadd = g + 0.2*cos(2*pi*t*180.327/512) + noise;

% freqz approach. Need to call with at least two frequencies
tic
h1 = freqz(gadd,1,2*pi*180.327/512*[1 1]);h1 = h1(1,1);
toc

Elapsed time is 0.042974 seconds.

timeit(@() freqz(gadd,1,2*pi*180.327/512*[1 1]) )

ans = 3.1574e-04

% sum approach
tic
h2 = sum(gadd.*exp(-1i*2*pi*t*180.327/512)); 
toc

Elapsed time is 0.004069 seconds.

timeit(@() sum(gadd.*exp(-1i*2*pi*t*180.327/512)))

ans = 2.2825e-05

sum is about an order of magnitude faster using either timeit or tic/toc. Curious as to why the timeit times are two orders of magnitude lower.

The results are pretty close to each other.

format long e
h1-h2

ans = 
     -6.750155989720952e-13 - 2.153832667772804e-12i

[xd]

Thanks.

I know the frequency since it is a stanard reference added , but it is dispersed in frequency domain since discrete calculation, so it is not at just right at digital coordinate of FFT.

But can not use 'interpolate' in frequency domain for FFT data.

So I don't know how to find it with my a little knowledge for FFT.

This give me a start point to use it for 2D --- that is my real problem.

In fact, you had answered my last post several days ago, but that time I do not know what is freq you suggested.

Wow, All you submit with comment disscuss mode but not answer mode, I unable select to accetpt the answer ！

[Star Strider]

‘Interpolating’ risks creeating data that did not previously exist.

If you zero-pad before calculating the Fourier transform, you increase the frequency resolution of the transform. Those will be the actual, calculated results at more frequencies, so no interpolation is necessary, and you do not create data where none previously existed.

Answer 1

[Paul]

Hi xd.

Here is one way to think about this problem.

Let c[n] = A*cos(w_0*n + phi) be the reference signal where -pi <= w_0 <= pi and -pi <= phi <= pi.

Let w[n] be a window function that is 0 for n < 0 and n >= N.

The Discrete Time Fourier Transform (DTFT) of g[n] = c[n]w[n] is

where W(w) is the DTFT of w[n].

Hence, if we know G(w) at some frequency we can, in principal, recover A and phi. However, that equation doesn't have a closed form solution for A and phi. But, if we choose a frequency where W(w - w0) is large and W(w + w0) is small, e.g., at w = w0, then we have an approximation

from which A and phi are readily obtained.

The DTFT of a rectangular window of length N is

(implemented this way this function can be a bit wonky for w very close to non-zero, integer multiples of 2*pi, but it's good enough for this exercise):

W = @(w,N) fillmissing((1-exp(-1j*w*N))./(1-exp(-1j*w)),'constant',N);
%W = @(w,N) sum(exp(-1j*w.*(0:N-1).'),1);  % straight implementation by definiton

For example

N   = 512;
w_0 = 2*pi*180.327/N;
A   = 1.8;
phi = 148.123*pi/180;
n   = 0:N-1;
g   = A*cos(w_0*n + phi);

Compute G(w0)

Gw_0 = sum(g.*exp(-1j*w_0*n));

Divide by the approximate contribution of the window

C = Gw_0/((W(0,N)/2));

Recover amplitude and phase and compare to truth

[abs(C),A]

ans = 1×2

1.8012 1.8000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[angle(C),phi]*180/pi

ans = 1×2

148.2411 148.1230

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

We see that for this case the effect of the window is basically negligible.

Let's plot the separate contributions of the window

w = (0:2^14-1)/2^14*2*pi;
figure
plot(w,abs(W(w-w_0,N)),w,abs(W(w+w_0,N))),grid
xline(w_0,'g')
xlabel('\omega (rad/sample)');ylabel('abs')
legend('W(\omega-\omega_0)','W(\omega+\omega_0)')

We see that W(w+w0) is very small at w = w_0.

What happens if we increase the frequency of the reference to something closer to pi rad/sample

w_0 = 2*pi*253.327/N;
g = A*cos(w_0*n + phi);
Gw_0 = sum(g.*exp(-1j*w_0*n));
C = Gw_0/((W(0,N)/2));
[abs(C),A]

ans = 1×2

1.7485 1.8000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[angle(C),phi]*180/pi

ans = 1×2

150.6970 148.1230

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

The estimates of A and phi are a bit off because the W(w+w0) contribution is not quite negligible at w = w_0 for the parameters of this problem (try instead setting w_0 = 2*pi*253/N as an experiment).

figure
plot(w,abs(W(w-w_0,N)),w,abs(W(w+w_0,N))),grid
xline(w_0,'g')
xlabel('\omega (rad/sample)');ylabel('abs')
legend('W(\omega-\omega_0)','W(\omega+\omega_0)')
xlim([-0.1,0.1]+w_0)

In this case, we can resort to solving for A and phi numerically

f1 = @(x) Gw_0 - (x(1))/2*(exp(1j*(x(2)))*W(w_0-w_0,N) + exp(-1j*(x(2)))*W(w_0+w_0,N));
f2 = @(x) [real(f1(x)),imag(f1(x))];
x = fsolve(f2,[abs(Gw_0)/(N/2),angle(Gw_0)]);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

[x(1),A]

ans = 1×2

1.8000 1.8000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[x(2),phi]*180/pi

ans = 1×2

148.1230 148.1230

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

If g[n] includes components in addition to the reference signal, e.g.,

g[n] = s[n]w[n] + c[n]w[n]

then by linearity the DTFT of g[n] is

where F(*) denotes the DTFT of the argument. If that first term on the right is small at a frequency where the second term on the right is large, such as w_0, then we can probably reasonably recover A and phi using the same procedure as above. If not (assuming the specific form of F(s[n]w[n]) is unknown), then there may be limitations as to how accurately A and phi can be estimated.

[xd]

Thanks，you are excellent expert !

I almost unstand you, suppose f1 is a liner function based on linearity of the DTFT, though this linearity maybe diffcult appear since Gw_0 is complicated in real data .

I am not good for frequency analysis, so I will learn to sure what is (1-exp(-1j*w*N))./(1-exp(-1j*w)

and what is exp(1j*(x(2)))*W(w_0-w_0,N) + exp(-1j*(x(2)))*W(w_0+w_0,N).

Of course，in real data, I will also need to check how about noise, error and calculate time 。