Interpolation of cross-correlation from zero padded idft of power spectrum shows odd behavior
13 次查看(过去 30 天)
显示 更早的评论
Hello,
I'm trying to interpolate the cross correlation between two singals by zero padding the ifft of the signals' PSD. Can anyone explain the odd spikiness of the interpolated correlation?
When the dft is evalauted directly via goertzel's method at non-integer indices, the same spikiness shows up.
Thanks for taking a look at this!
Signals being correlated
d1 = designfilt("lowpassiir",FilterOrder=12, ...
HalfPowerFrequency=0.02,DesignMethod="butter");
N = 2^10;
y1 = rand(N,1);
y1 = filtfilt(d1,y1)';
N2 = 2^2;
shiftImposed = 1*N2;
y2 = interp1((1:N),y1,(1+shiftImposed:N+shiftImposed));
y3 = imresize(y1,1/N2);
y4 = imresize(y2,1/N2);
y5 = y3(length(y3)/8:length(y3)-length(y3)/8+1);
y6 = y4(length(y3)/8:length(y3)-length(y3)/8+1);
figure(54);plot(y5);hold on;plot(y6)
Correlation without interpolation
g1 = y5-mean(y5);
g2 = y6-mean(y6);
G1 = fft(g1);
G2 = fft(g2);
R = G1.*conj(G2);
N = length(R);
p=0;
N2 = 2^p*N;
x = (1:N2);
R2 = padarray(R,[0 (N2-N)/2] ,0,'both');
rr2 = fftshift(ifft(R2,'symmetric'));
figure(43);plot(abs(rr2))
Correlation with 2N interpolation via zero padding ifft of cross power spectrum
g1 = y5-mean(y5);
g2 = y6-mean(y6);
G1 = fft(g1);
G2 = fft(g2);
R = G1.*conj(G2);
N = length(R);
p=1;
N2 = 2^p*N;
x = (1:N2);
R2 = padarray(R,[0 (N2-N)/2] ,0,'both');
rr2 = fftshift(ifft(R2,'symmetric'));
figure(43);plot(abs(rr2))
Goertzel's Method to calculate dft bins of arbitrary location on unit circle shows the same thing as figure 3
x = (1:.5:N);
DD = fftshift(goertzel(R,x));
figure(457);plot(abs(DD))
0 个评论
采纳的回答
Paul
2022-12-1
Hi matt,
I think the issue is with how the padding is being done. See code below for that and other questions. Full disclosure: I'm only looking at how the fft/ifft is being used; I'm not familiar with this type of problem.
rng(100);
d1 = designfilt("lowpassiir",FilterOrder=12, ...
HalfPowerFrequency=0.02,DesignMethod="butter");
N = 2^10;
y1 = rand(N,1);
y1 = filtfilt(d1,y1)';
N2 = 2^2;
shiftImposed = 1*N2;
y2 = interp1((1:N),y1,(1+shiftImposed:N+shiftImposed));
y3 = imresize(y1,1/N2);
y4 = imresize(y2,1/N2);
y5 = y3(length(y3)/8:length(y3)-length(y3)/8+1);
y6 = y4(length(y3)/8:length(y3)-length(y3)/8+1);
figure(54);plot(y5);hold on;plot(y6)
Correlation without interpolation
g1 = y5-mean(y5);
g2 = y6-mean(y6);
G1 = fft(g1);
G2 = fft(g2);
R = G1.*conj(G2);
N = length(R);
p=0;
N2 = 2^p*N;
x = (1:N2);
% no padding here
N2 - N
R2 = padarray(R,[0 (N2-N)/2] ,0,'both');
Why is fftshift applied here?
rr2 = fftshift(ifft(R2,'symmetric'));
figure(43);plot(abs(rr2))
Correlation with 2N interpolation via zero padding ifft of cross power spectrum
g1 = y5-mean(y5);
g2 = y6-mean(y6);
G1 = fft(g1);
G2 = fft(g2);
R = G1.*conj(G2);
N = length(R);
p=1;
N2 = 2^p*N;
x = (1:N2);
% pad with this many 0's in both directions
N2-N
I think the issue is here. Need to fftshift R, and then pad the left and the right for negative and positive frequencies, and then ifftshift back for use with ifft in the next step.
%R2 = padarray(R,[0 (N2-N)/2] ,0,'both');
R2 = ifftshift(padarray(fftshift(R),[0 (N2-N)/2] ,0,'both'));
Why is fftshift applied here?
rr2 = fftshift(ifft(R2,'symmetric'));
figure(43);plot(abs(rr2))
0 个评论
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!