Remove peak at 0 hz of fft

조회 수: 130 (최근 30일)
Emmy
Emmy 2020년 7월 6일
댓글: Image Analyst 2020년 7월 10일
Hi everybody,
I want to calculate the fft of my data, but I get a large peak at 0 Hz. I have read all the answers to the previous questions. Mostly recommended using x=x-mean(x) or x=detrend(x). I have tried that but nothing works.
Another suggestion was to look at the diff(x) and diff(diff(x)). I tried that and both give a noisy straight line.
I dont now what else to do so can somebody help me?
Thanks in advance!
x=load('data.txt');
% diffx=diff(x);
% diffdiffx=diff(diff(x));
xm=x-mean(x);
xd=detrend(x);
y=detrend(xm,'constant');
N=length(x);
fs=1000;
Ts=1/fs;
f=(0:1/N:1/2-1/N)*fs;
% figure
% plot(x)
% hold on
% plot(diffx)
% plot(diffdiffx)
% hold off
% legend('x','diffx','diffdiffx')
xhat=fft(x,N);
xhat(1)=0;
figure
plot(f,abs(xhat(1:N/2)));
title('Normal data')
xhat=fft(xm,N);
xhat(1)=0;
figure
plot(f,abs(xhat(1:N/2)));
title('Removed mean data')
xhat=fft(xd,N);
xhat(1)=0;
figure
plot(f,abs(xhat(1:N/2)));
title('Detrend data')
xhat=fft(y,N);
xhat(1)=0;
figure
plot(f,abs(xhat(1:N/2)));
title('Removed mean and detrend data')
[S,fm]=periodogram(x,[],N,fs,'power');
figure
plot(fm,10*log10(S))
title('Periodogram')
[S,fm]=periodogram(xm,[],N,fs,'power');
figure
plot(fm,10*log10(S))
title('Periodogram removed mean')
[S,fm]=periodogram(xd,[],N,fs,'power');
figure
plot(fm,10*log10(S))
title('Periodogram detrend data')
[S,fm]=periodogram(y,[],N,fs,'power');
figure
plot(fm,10*log10(S))
title('Periodogram removed mean and detrend data')
[Sw,fw]=pwelch(x,ones(N,1),0,N,fs,'power');
figure
plot(fw,Sw)
title('Welch method')
[Sw,fw]=pwelch(xm,ones(N,1),0,N,fs,'power');
figure
plot(fw,Sw)
title('Welch method removed mean')
[Sw,fw]=pwelch(xd,ones(N,1),0,N,fs,'power');
figure
plot(fw,Sw)
title('Welch method detrend data')
[Sw,fw]=pwelch(y,ones(N,1),0,N,fs,'power');
figure
plot(fw,Sw)
title('Welch method removed mean and detrend data')
  댓글 수: 2
Mehmed Saad
Mehmed Saad 2020년 7월 6일
apply a hpf
Emmy
Emmy 2020년 7월 6일
Applying a highpass filter did not solve the problem. There still is a peak from the point of the high pass filter. And it doesn't matter at which point I take the high pass filter. Even when I look in the normal data, where the peak ends, there still is a peak after the high pass filter.

댓글을 달려면 로그인하십시오.

채택된 답변

David Goodmanson
David Goodmanson 2020년 7월 9일
Hi Emmy
It's all in how you view things. Literally. Here is what happens when using x - mean(x).
fs = 1000;
N = 30000;
figure(1)
plot(x)
xm = x - mean(x);
y = fftshift(abs(fft(xm))); % shift zero frequency to the center
f = (-N/2:N/2-1)*fs/N; % make f array match
y(15001) = 1; % adjust 0 Hz value in fft
figure(2)
semilogy(f,y)
figure(3)
semilogy(f,y)
xlim([-50 50])
Both x and y have 30000 elements. Ordinarily, y(1) would be the element corresponding to zero frequency. For plotting purposes, fftshift puts that element at the center of the array, element 15001.
Since the y values have a large range, it makes sense to look at abs(y) on a semilog plot. Note that since y is the fft of a real function, abs(y) is symmetric about f==0.
After subtracting the mean, the DC value should be zero. Because of 'numerical error' it is small but nonzero. For plot purposes I arbitrarily set that value to 1 so as to not have a large negative spike in the semilog plot.
In figure 2, you can see some lower frequency stuff rising above the noise. In the zoomed-in figure 3. there is something interesting happening between about 17-20 Hz. If you zoom in a little bit in figure 3 you will see no peak at f=0, with the largest values at one freq point on either side, falling off from there.
  댓글 수: 1
Image Analyst
Image Analyst 2020년 7월 10일
Adding plots for visualization purposes:
x=load('data.txt');
fs = 1000;
N = 30000;
subplot(3, 1, 1);
plot(x)
grid on;
title('x', 'FontSize', 20);
xm = x - mean(x);
y = fftshift(abs(fft(xm))); % shift zero frequency to the center
f = (-N/2:N/2-1)*fs/N; % make f array match
y(15001) = 1; % adjust 0 Hz value in fft
subplot(3, 1, 2);
semilogy(f,y)
grid on;
title('FT', 'FontSize', 20);
subplot(3, 1, 3);
semilogy(f,y)
xlim([-50 50])
grid on;
title('Zoomed FT', 'FontSize', 20);

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Mehmed Saad
Mehmed Saad 2020년 7월 6일

FFT of X

FFT of X High pass filtered

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by