What is the cause of strange "spike" artifacts after IFFT?

Question

Darcy 2025-7-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178683-what-is-the-cause-of-strange-spike-artifacts-after-ifft

编辑： Star Strider 2025-7-25

在 MATLAB Online 中打开

I am convolving two signals by multiplying them in the frequency domain. The first signal,

, is in the time domain, and I convert it to the frequency domain using an fft to get

. Prior to fft, I precondition the signal by detrending, applying a tukey window with 10% taper, and padding with zeros. The second signal is already in the frequency domain,

, so the only thing I need to do is interpolate B onto the frequencies defined in A.

I then multiply them to get

. To convert

to the time domain, I use an ifft to get

and then remove the padding.

The resulting

generally looks okay, except it has these very suspicious spikes in the time domain. See the attached files which show the full time series and then a couple zoom views of one of the spikes.

Any ideas what could be causing this? See example code below and example dataset for

and

.

clearvars; close all
pad = 100000;
load('at.mat');
load('Bf.mat')
%Bf(1:10) = [];
%freq(1:10) = [];
L = length(at); % number of time samples
t = 0:1/fs:(1/fs)*L-1/fs; %time vector
% ------------------------Preconditioning a(t)---------------------------
%Detrend a(t)
at = detrend(at);
%Apply tukey window
tukey = tukeywin(length(at),0.2);
at = tukey.*at;
%Add zero padding
at = cat(1,ones(pad,1).*at(1),at,ones(pad,1).*at(end));%-nanmean(xt);
Lpad = length(at);
%--------------------Convert a(t) to A(f)--------------------------------
%Compute a(t) in frequency domain A(f)
Af = fft(at); 
df = fs/Lpad; %frequency resolution
fAxis = (0:df:(fs-df)) - (fs-mod(Lpad,2)*df)/2; %full frequency axis including negative freqs
fAxis(abs(fAxis)<df/1000000) = 0;
%Get positive frequencies only, careful to consider odd and even lengths
if mod(Lpad,2)==0
    f = fAxis((Lpad/2)+1:end);
else
    f = fAxis((Lpad+1)/2:end);
end
%-------------------------------------------------------------------------
%-------------------Interpolating B(f)----------------------------------
%B(f) is imaginary, so interpolate real and imaginary separately
realB = interp1(freq,real(Bf),abs(fAxis(fAxis<0)),'linear')';
imagB = interp1(freq,imag(Bf),abs(fAxis(fAxis<0)),'linear')';
%Reconstruct interpolated complex B(f)
Bq = realB+1i*imagB;
%Ensure conjugate symmetry of B(f) prior to computing C(f)
if mod(length(fAxis),2)==0
    Bconj = fftshift(cat(1,conj(Bq),0,flip(Bq(2:end))));   
else
    Bconj = fftshift(cat(1,conj(Bq),0,flip(Bq)));  
end
%-----------------Compute C(f)---------------------------------------
Cf = Af.*Bconj;
%If the query frequencies extend beyond the frequency range of B(f), then
%the interpolate sets them to NaN. These are set to zero, rather than
%removing them because removing them alters the number of points in the fft
%and ifft
Cf(isnan(Cf))=0;
ct = ifft(Cf,'symmetric');
%Remove padding
ct(1:pad) = [];
ct(end-pad+1:end) = [];
figure
plot(t,ct,'-k');

2 个评论
显示无隐藏无

Matt J 2025-7-23

Any ideas what could be causing this?

No, we don't know what code you ran.

Darcy 2025-7-23

编辑：Darcy 2025-7-23

@Matt J Updated to include a code example and data sample

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2025-7-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178683-what-is-the-cause-of-strange-spike-artifacts-after-ifft#answer_1568368

The signals are approximately the same lengths in time, however 'Zoom1' has a significantly higher sampling frequency than 'Zoom2'. It appears from 'Full_TimeSeries_c' that you interpolated 'Zoom2' to match 'Zoom1'. This has the effect of creating data where no data previously existed, so the interpolation function is simply guessing what those data should be.

It would likely be better (at least more justifiable) to interpolate 'Zoom1' to match 'Zoom2'. That still risks creating data, however those data are well-described in the original signal, so those interpolated data at least have some basis in reality. You are essentially just using fewer data from 'Zoom1', rather than creating 'Zoom2' data that do not actually exist.

At least try that to see if the result improves.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Star Strider 2025-7-24

在 MATLAB Online 中打开

The 'spikes' appear to be in 'Af' in the frequency domain, and they are likely propagated throughout your code. With respect to 'interpolation', interpolating a 38-sample signal out to 2,273,600 samples is a bit of a stretch, to say the least. You have no idea what 'Bf' is actually doing at those sample points, for example it could be oscillating wildly and you sould simply never catch that behaviour.

I suggest that you reverse this approach and interpolate 'at' to the size of 'Bf', in whatever domain you choose. I cannot determine if that will improve the result, however it definitely makes more sense (at least to me).

clearvars; close all

pad = 100000;

load('at.mat');

load('Bf.mat')

whos

Name Size Bytes Class Attributes Bf 38x1 608 double complex at 2073600x1 16588800 double freq 38x1 304 double fs 1x1 8 double pad 1x1 8 double

fs

fs = 8

pad

pad = 100000

figure

plot(at)

grid

title('at')

figure

plot(freq, real(Bf), '.-', DisplayName='Re(B_f)')

hold on

plot(freq, imag(Bf), '.-', DisplayName='Im(B_f)')

plot(freq, abs(Bf), '.-', DisplayName='|B_f|')

hold off

xlabel('Frequency')

ylabel('Amplitude')

legend(Location='best')

title('B_f')

%Bf(1:10) = [];

%freq(1:10) = [];

L = length(at); % number of time samples

t = 0:1/fs:(1/fs)*L-1/fs; %time vector

% ------------------------Preconditioning a(t)---------------------------

%Detrend a(t)

at = detrend(at);

%Apply tukey window

tukey = tukeywin(length(at),0.2);

at = tukey.*at;

%Add zero padding

at = cat(1,ones(pad,1).*at(1),at,ones(pad,1).*at(end));%-nanmean(xt);

Lpad = length(at);

%--------------------Convert a(t) to A(f)--------------------------------

%Compute a(t) in frequency domain A(f)

Af = fft(at);

df = fs/Lpad; %frequency resolution

fAxis = (0:df:(fs-df)) - (fs-mod(Lpad,2)*df)/2; %full frequency axis including negative freqs

fAxis(abs(fAxis)<df/1000000) = 0;

%Get positive frequencies only, careful to consider odd and even lengths

if mod(Lpad,2)==0

f = fAxis((Lpad/2)+1:end);

else

f = fAxis((Lpad+1)/2:end);

end

%-------------------------------------------------------------------------

%-------------------Interpolating B(f)----------------------------------

%B(f) is imaginary, so interpolate real and imaginary separately

realB = interp1(freq,real(Bf),abs(fAxis(fAxis<0)),'linear')';

imagB = interp1(freq,imag(Bf),abs(fAxis(fAxis<0)),'linear')';

% [realB,fr] = resample(real(Bf),freq,1/mean(diff(abs(fAxis(fAxis>0)))),'linear');

% [imagB,fi] = resample(imag(Bf),freq,1/mean(diff(abs(fAxis(fAxis>0)))),'linear');

%Reconstruct interpolated complex B(f)

Bq = realB+1i*imagB;

%Ensure conjugate symmetry of B(f) prior to computing C(f)

if mod(length(fAxis),2)==0

Bconj = fftshift(cat(1,conj(Bq),0,flip(Bq(2:end))));

else

Bconj = fftshift(cat(1,conj(Bq),0,flip(Bq)));

end

SAf = size(Af)

SAf = 1×2

2273600 1

SBconj = size(Bconj)

SBconj = 1×2

2273600 1

%-----------------Compute C(f)---------------------------------------

Cf = Af.*Bconj;

figure

plot(real(Af), DisplayName='Re(Af)')

hold on

plot(imag(Af), DisplayName='Im(Af)')

plot(real(Bconj),':', DisplayName='Re(B_{conj})')

plot(imag(Bconj),':', DisplayName='Im(B_{conj})')

hold off

ylim([-1 1]*1E3)

legend(Location='best')

title('A_f & B_{conj}')

figure

plot(real(Af), DisplayName='Re(Af)')

hold on

plot(imag(Af), DisplayName='Im(Af)')

plot(real(Bconj),':', DisplayName='Re(B_{conj})')

plot(imag(Bconj),':', DisplayName='Im(B_{conj})')

hold off

ylim([-1 1]*1E3)

xlim([0.27 0.45]*1E6)

legend(Location='best')

title('A_f & B_{conj} (Zoomed)')

%If the query frequencies extend beyond the frequency range of B(f), then

%the interpolate sets them to NaN. These are set to zero, rather than

%removing them because removing them alters the number of points in the fft

%and ifft

Cf(isnan(Cf))=0;

ct = ifft(Cf,'symmetric');

%Remove padding

ct(1:pad) = [];

ct(end-pad+1:end) = [];

figure

plot(t,ct,'-k');

.

Star Strider 2025-7-24

编辑：Star Strider 2025-7-25

在 MATLAB Online 中打开

If you have the Signal Processing Toolbox, a better way to do this is to use the invfreqz function to create a digital filter from 'Bf' and then use filtfilt to filter 'at'.

There are significant problems with your data. To begin with, 'Bf' does not have regularly-spaced frequencies. Essentially everything in signal processing requires equally-spaced independent variable values. In order for the filter to work correctly, 'at' needs to have the same sampling frequency as 'Bf' and the original 'Bf' time-domain signal is defined at being sampled as two times the Nyquist frequency. The only reference we have for that is the highest frequency for 'Bf' (1.7813) so the sampling frequency is 3.5625. (The units aren't stated, so I just assume they're compatible.) In order for the filter to work, 'at' needs to be resampled at that frequency, and that was done with the resample function. (I use that twice here, once for 'Bf' in the frequency domain, and once for 'at' in the time domain. I just appended 'r' to the respective variable names to create the new vectors.)

Try something like this --

load('at.mat')

load('Bf.mat')

L = length(at); % number of time samples

t = 0:1/fs:(1/fs)*L-1/fs; %time vector

Fs = 1/mean(diff(t))

Fs = 8

disp(freq)

1.7813 1.3437 1.0000 0.7500 0.5625 0.4062 0.3125 0.2344 0.1719 0.1328 0.1016 0.0781 0.0586 0.0430 0.0320 0.0242 0.0180 0.0133 0.0102 0.0076 0.0057 0.0043 0.0033 0.0025 0.0018 0.0013 0.0010 0.0007 0.0005 0.0004 0.0003 0.0002 0.0002 0.0001 0.0001 0.0001 0.0001 0.0000

% dfreq = diff(freq)

Fsfreq = mean(diff(freq))

Fsfreq = -0.0481

Fss = -1/Fsfreq % Mean Frequency Sampling-frequency for 'Bf'

Fss = 20.7723

Bffs = 2*max(freq) % 'Bf' Sampling Frequency (2x Nyquist Frequency)

Bffs = 3.5625

[atr,tr] = resample(at, t, Bffs); % Resample 'at' With 'Bf' Sampling Frequency

[Bfr,freqr] = resample(Bf, freq, Fss); % Resample 'Bf' To Regularly-Spaced Frequencies

% [Bfr,freqr] = resample(Bf, freq, fs)

figure

stem(freq,'.-')

grid

% figure

% stem(dfreq,'.-')

% grid

figure

plot(flip(freq), abs(flip(Bf)),'.-')

grid

title('Original ''B_f''')

figure

plot(freqr, abs(Bfr), '.-')

grid

title('Resampled ''B_f''')

iBf = imag(Bf);

[pks,plocs] = findpeaks(iBf);

ixz = find(diff(sign(iBf)));

for k = 1:numel(ixz)

idxrng = max(1,ixz(k)-1) : min(numel(freq),ixz(k)+1);

fv(k) = interp1(iBf(idxrng), freq(idxrng), 0);

% iBfv(k) = interp1(freq(idxrng), iBf(idxrng), fv(k))

end

[fvu,ia] = uniquetol(fv, 0.1);

figure

plot(freq, imag(Bf),'.-')

hold on

plot(freq(plocs), pks, '^r')

plot(fvu, zeros(size(fvu)), 'og')

hold off

grid

xlabel('Frequency')

ylabel('Magnitude')

title('Original Im(''B_f'')')

yline(0, ':k', LineWidth= 1.5)

text(0.1, 0.02, ["Poles = 3","Zeros = 3","Zero at the origin","Possible Pole at +Inf"])

iBfr = imag(Bfr);

[pksr,plocsr] = findpeaks(iBfr);

yl = 0.0001;

ixz = find(diff(sign(iBfr-yl)));

for k = 1:numel(ixz)

idxrng = max(1,ixz(k)-1) : min(numel(freqr),ixz(k)+1);

fvr(k) = interp1(iBfr(idxrng), freqr(idxrng), yl);

% iBfvr(k) = interp1(freqr(idxrng), iBfr(idxrng), fvr(k))

end

fvr

fvr = 1×6

0.0031 0.2803 0.3078 0.5480 0.5480 1.4021

[fvru,ia] = uniquetol(fvr, 0.1)

fvru = 1×4

0.0031 0.2803 0.5480 1.4021

●

ia = 4×1

1 2 4 6

% fvru = fvr;

figure

plot(freqr, imag(Bfr),'.-')

hold on

plot(freqr(plocsr), pksr, '^r')

plot(fvru, zeros(size(fvru))+yl, 'og')

hold off

grid

xlabel('Frequency')

ylabel('Magnitude')

title('Resampled Im(''B_f'')')

yline(yl, ':k', LineWidth= 1.5)

text(0.1, 0.02, ["Poles = 3","Zeros = 4","Zero at the origin","Possible Pole at +Inf"])

figure

plot(freqr, abs(Bfr),'.-')

grid

xlabel('Frequency')

ylabel('Magnitude')

title('Resampled ''B_f'' (Absolute Value)')

[b,a] = invfreqz(Bf, freq, 4, 4, [], 100);

figure

freqz(b, a, 2^16, Bffs)

% sgtitle('Filter Response')

set(subplot(2,1,1), 'Ylim',[min(ylim) 0], 'XScale','log')

set(subplot(2,1,2), 'XScale','log')

[br,ar] = invfreqz(Bfr, freqr, 4, 4, [], 100);

figure

freqz(br, ar, 2^16, Bffs)

% sgtitle('Filter Response')

set(subplot(2,1,1), 'Ylim',[min(ylim) 0], 'XScale','log')

set(subplot(2,1,2), 'XScale','log')

atr_filt = filtfilt(br, ar, atr);

[FTs1,Fv] = FFT1(atr,tr);

[FTsf1,Fvf] = FFT1(atr_filt,tr);

figure

tiledlayout(2,1)

nexttile

semilogx(Fv, abs(FTs1))

grid

xlim('tight')

xlabel('Frequency')

ylabel('Magnitude')

title('Resampled ''at''')

nexttile

semilogx(Fvf, abs(FTsf1))

grid

xlim('tight')

xlabel('Frequency')

ylabel('Magnitude')

title('Resampled Filtered ''at''')

figure

tiledlayout(2,1)

nexttile

plot(tr, atr)

grid

xlim('tight')

ylim([1.375 1.55]*1E+4)

title('Resampled ''at'' (''at_r'')')

xlabel('Time')

ylabel('at_r')

nexttile

plot(tr, atr_filt)

grid

xlim('tight')

ylim([0.03 0.034])

xlabel('Time')

ylabel('at_r (filtered)')

function [FTs1,Fv] = FFT1(s,t)

% One-Sided Numerical Fourier Transform

% Arguments:

% s: Signal Vector Or Matrix

% t: Associated Time Vector

t = t(:);

L = numel(t);

if size(s,2) == L

s = s.';

end

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));

Fv = Fs*(0:(NFFT/2))/NFFT;

% Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

Fv = Fv(:);

FTs1 = FTs(Iv,:);

end

The derived filter (with all the necessary corrections required by signal processing conventions) actually has very little effect on 'atr' because it is essentially negligible in the frequencies that appear to be important to describing 'atr'.

The 'spikes' were an artifact of the original signal processing approach. They are not present when all the necessary manipulations of your signals that signal processing requires are met.

EDIT -- (24 Jul 2025 at 16:40)

Added explanations and corrected code to conform to accepted (and essentially required) signal processing conventions.

EDIT -- (25 Jul 2025 at 00:30)

Added (restored) the '[b,a] = invfreqz(Bf, freq, 4, 4, [], 100)' code to be complete. It likely does not make any difference which filter you use ('[b,a]' or '[br,ar]') since the passband is sufficiently high that it will only accentuate the 1 (Hz?) noise, and will not affect the 'atr' signal at all, except to attenuate it.

.

请先登录，再进行评论。

What is the cause of strange "spike" artifacts after IFFT?

2 个评论
显示无隐藏无

回答（1 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

What is the cause of strange "spike" artifacts after IFFT?

2 个评论 显示 无隐藏 无

回答（1 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论