How to remove DC component in FFT?

I succesfully plotted my FFT with MATLAB discussion help. Now I could not remove the DC component at 0Hz. Which shows me a very high amplitude. Can any one suggest me an idea?

data1 = xlsread('Reading 1.xlsx') ; %Loading Sensor data from Excel file
t = data1 (1:512,2);                %Selecting Time vector
s = data1 (1:512,3);                %Selecting Z axis vibrations
L = numel(t);                       %Signal length
Ts = mean(diff(t));                 %Sampling interval
Fs = 1/Ts;                          %Sampling frequency
Fn = Fs/2;                          %Nyquist frequency
FTs = fft(s)/L;                     %Fast fourier transform (s- data)
Fv = linspace(0,1, fix(L/2)+1)*Fn;  %Frequency vector
Iv = 1:numel(Fv);                   %Index vector

subplot(2, 1, 1);                   %plotting top pane
plot(t,s);                          %Acceleration vs time
set(gca,'xlim',[1 50]);             %Scale to fit
grid;                               %Grids on
title ('Acceleration vs time');     
xlabel('time(s)');
ylabel('Acceleration');

subplot(2, 1, 2);                   %Plotting bottom pane
plot(Fv, abs(FTs(Iv))*2,'red');     %FFT - Amplitude vs Frequency
grid
title ('Fast fourier transform'); 
xlabel('Frequency (Hz)');
ylabel ('Amplitude (m)');

Answer 1

[Image Analyst]

In the spatial domain, before fft, you can subtract the mean

Iv = Iv - mean(Iv);

In the frequency domain, you can zero out the DC component by setting it to zero

ft = fft(Iv);
ft(1) = 0;

[Wakeel Mohammed]

First of all thank you so much for taking your time on my question.

Lemme try this and give a positive feedback!

[Nidhi Singh]

Have u tried this code? I am not getting zero out the dc offset by only making it component zero. Is there any other method to remove high amplitude at zero frequency.

[Image Analyst]

@Nidhi Singh, it should work, or at least be very close.

ft = fft(timeDomainSignal); % Tranform the original signal.
ft(1) = 0; % Eliminate DC component so there will be no offset to the signal in the time domain.
timeDomainSignal = ifft(ft); % Inverse transform back to time domain.
% Compute the mean.  It should now be close to zero.
meanValue = mean(real(timeDomainSignal))

The mean value should be pretty close to zero.

[Cédric Cannard]

Thanks for this @Image Analyst!

[Sateesh Kandukuri]

Dear @Image Analyst, I have tried all your methods to remove the dc component for my data, but nothing worked.

I have tried in all ways from different posts, but ended up with the same results.

clear;

A = readmatrix('table.txt');

ts=1e-12;

Mz = A(:,4);

Mz = Mz - mean(Mz);

N = 2^(nextpow2(length(Mz)));

freq = fft(Mz,N);

freq(1)=0;

newMz = ifft(freq);

newMz = newMz-mean(newMz);

freq2 = abs(fftshift(fft(newMz,N)));

F = ((-N/2:N/2-1)/N)*(1/ts);

plot(F(N/2+1:end)/1e9,freq2(N/2+1:end),'b','Linewidth',1);

Could you please check where I am doing wrong. I attached the table for your reference.

[Image Analyst]

It is removed. In fact you triple removed it. You subtracted the mean (removal #1) then you fft'd it and set the first element to 0 (removal #2) and then you ifft'd is and again subtracted the mean (removal #3). You only need to do one of those, not all three.

However after your round trip to the spatial domain you take the absolute value of the signal. So instead of the signal going positive and negative, which will give the mean of 0 (no DC component), you flipped all the negative values to positive. Now there is no way the mean will be zero anymore unless the entire signal is all zero because you have tons of positive, non-zero values. Don't call abs().

[Sateesh Kandukuri]

Dear @Image Analyst, Thank you so much for your quick response. Actually, in the programme, I kept all the ways I approached removing the dc component.

I have tried without the abs() function, but the high amplitude peak at zero frequency was not removed.

clear;

A = readmatrix('table.txt');

ts=1e-12;

Mz = A(:,4);

Mz = Mz - mean(Mz);

N = 2^(nextpow2(length(Mz)));

freq = abs(fftshift(fft(Mz,N)));

% freq = freq/max(freq);

freq(1)=0;

F = ((-N/2:N/2-1)/N)*(1/ts);

plot(F(N/2+1:end)/1e9,freq(N/2+1:end),'b','Linewidth',1);

xlabel('Frequency (GHz)')

ylabel('Spectrum (a.u.)')

xlim([0 20]);

With and without the abs() function I got the following results.

[Image Analyst]

Glad it is working now. The DC component is now removed, but as you can see there is power in frequencies near pure zero. You may want to erase a longer section at the beginning and end of your unshifted spectrum array to get rid of those. Like

freq(1:30) = 0;
freq(end-30 : end) = 0;

Or you can do it in the middle of the array if you're dealing with the shifted array.

[AMOS]

ft = fft(Iv);

ft(1) = 0;

Hi, how can the above code be implemented, in terms of ft as there is no ft in the initial code?

[AMOS]

At what point is ft implemented?

[Image Analyst]

@AMOS, at the point where you run this line of code:

ft = fft(Iv);

[Ajith Kumar]

Iv_new = (Iv / mean(Iv)) - 1

This removes the DC component by normalizing the signal with its mean and then centering it around zero.

Is this another possible solution?

Answer 2

[Sateesh Kandukuri]

Dear @Image Analyst, the behaviour of my system is asymmetrical to the excitation fields.

Is it possible to modify this behaviour from asymmetrical to symmetrical? And then performing FFT may resolve my issue.

[Image Analyst]

What would you hope the final signal would look like? And, perhaps more importantly, WHY do you need the signal to look like that?

[Sateesh Kandukuri]

@Image Analyst, The behaviour of my system without the excitation field in the relaxation process is like below

And the corrsponding fft is

But when I excite the system, the response follows the path of relaxation. That you can see clearly in the previous attached image. When I do the fft calculation of the excited system, the peak close to zero frequency is dominating clearly due to this issue. I can't make the initial and final points of fft to zero due to other peaks appearing close to those points.

[Image Analyst]

Why not juse use movmean() on the original signal with a window width that is the length of one of the small wave periods?

[Sateesh Kandukuri]

@Image Analyst, Could you please send the snippet of the code?

[Image Analyst]

Not sure how to get your signal. Attach x and y in a .mat file. It would be something like

windowWidth = 41; % Whatever
ySmooth = movmean(y, windowWidth);

[Sateesh Kandukuri]

@Image Analyst, I failed to use movmean() function correctly. I attached the input file for verification.

[Image Analyst]

Try this:

s = load('t_my.mat')
x = s.A(:, 1);
y = s.A(:, 2);
plot(x, y, 'b-')
grid on;
[peakValues, indexesOfPeaks] = findpeaks(y)
hold on;
plot(x(indexesOfPeaks), peakValues, 'rv')
windowWidth = 2 * mean(diff(indexesOfPeaks))
ySmooth = movmean(y, windowWidth);
plot(x, ySmooth, 'r-', 'LineWidth', 2)

[Sateesh Kandukuri]

@Image Analyst, It's working. You are a true lifesaver! Million thanks Sir.

[Sateesh Kandukuri]

Dear @Image Analyst, this logic is working for smoothly varying magnetization components(Mx and My). But in the case of Mz componet I got the following result.

I attached the input file for you to look over.

[Image Analyst]

I had the window width be several wavelents long. How many indexes are between each of your peaks? Try having the window width be like 3 or 4 times that long.

[Sateesh Kandukuri]

Dear @Image Analyst,

Using your suggestion, I used movmean() in the calculation of fft as

A = readmatrix('table.txt');

ts=1e-12;

My = A(:,3);

[peakValues, indexesOfPeaks] = findpeaks(My);

windowWidth = 2 * mean(diff(indexesOfPeaks));

MySmooth = movmean(My, windowWidth);

My = My - MySmooth;

N = 2^(nextpow2(length(My)));

freq = fft(My,N);

freq2 = abs(fftshift(freq));

freq3 = freq2/max(freq2);

I got the following result for My component

Is this the right way to use movmean() function?

I've tried to understand the working of movmean() function using some arrays, but I still need clarification. Can you briefly explain with an example?