Remove peak at 0 hz of fft
100 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Emmy
el 6 de Jul. de 2020
Comentada: Image Analyst
el 10 de Jul. de 2020
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 comentarios
Respuesta aceptada
David Goodmanson
el 9 de Jul. de 2020
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 comentario
Image Analyst
el 10 de Jul. de 2020
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);
Más respuestas (1)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!