Find the problem with my FFT code

Hello

I'm using FFT code to compare the amplitudes of two signals, but the answer is different from my idea.

I measured the pressure signal 1000 times a second.

And one is unchanged, the other is increasing.

I divided these two signals by 8 seconds each and PLOT through FFT.

If I find the amplitude with the code below, the value of the unchanging signal is bigger.

But in my opinion, the amplitude of the magnitude of the increasing signal should be bigger den unchanged signal.

Is my code wrong?

thank you for reading :)

Fs = 1000; % Sampling frequency
T = 1/Fs;  % Sample time
L = 8000;  % Length of signal
t = (0:L-1)*T;                 % Time vector
x=TL26;                        %signal
y=table2array(x);
n=4; %% order of Butterworth filter
Wn=[130 150]; % 130 ~ 150Hz signal
Fn=Fs/2; % Nyquist
ftype='bandpass';
[b,a]=butter(n, Wn/Fn, ftype);%%% butter(): 버터워스 returns an order 2*n digital bandpass filter if
%%%fbn is a two-element vector
e=filtfilt(b,a,y);
N = length(e);
NFFT = 2^nextpow2(N); % Next power of 2 from length of y
%Y = fft(y,NFFT)/NFFT;
Y = fft(e,NFFT);
Ya = abs(Y)/NFFT; % correctly normalised amplitude
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure(99)
plot(f,2*abs(Ya(1:NFFT/2+1)));
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

[dpb]

Save us some effort and show the results you obtained. Other than not correcting for DC component, I don't see anything obvious in the code...looking at the two time traces, however, it appears to my eye there's more regular frequency component in the SS case than in the transient one...it may be that's what the FFT is revealing.

Remember amplitude in the frequency space is on a per-frequency basis; the energy in the time series contains all frequency components so the amplitude encompasses all frequencies that contain any energy--only the energy in each frequency which is a fraction of the total is in each energy bin. So, if there is more of a given frequency in one signal than another even though the overall amplitude may be larger, it may be spread across more/different frequencies in frequency space.

[dongwoo kim]

Thank you for your answer.

So, what more should I do? I know that these two signals need to be FFT simply.

Can you give me examples or hints?

Thanks

[dpb]

Don't know there's anything you can do. The frequency content is what it is.

You didn't attach the two FFT plots and I don't feel like taking time to actually download and run...

[dongwoo kim]

Hi I'll show you two plots

Time plotted values between 100 and 108 seconds.

The large amplitude is the value where the pressure does not change.

The graph I first presented shows that the pressure value between 100 and 108 seconds is small but the plotted value is large.

[dpb]

Dunno how you jump to that conclusion -- the two peaks are at roughly 0.0009 signal/Hz and differ by <5%. That's well less than sqrt(Ypk) the standard error of estimate of single estimate. By comparison, the RMS of the time signal is 0.6--roughly 600X the peak amplitude. Eve

The RMS level of the time signals is also essentially the same...there's a little more lower frequency beats in the R signal than T and vice versa for some higher frequencies showing up for whatever reason.

But, I certainly see nothing suspicious at all in those figures...

[dongwoo kim]

I compared the other signal files from 100 to 108 seconds.

The obviously unchanged signal is the y value.

[dpb]

And? What do you expect? Looks perfectly reasonable...

[dongwoo kim]

I knew if I showed a new graph, I could see the wrong part of my code. I agree with you that the margin of error is too narrow.

[dpb]

What wrong part of code is that?

[dongwoo kim]

As I said before, I think the fft size of the growing signal should be large. But small. This is a little difference, but it's not like I thought the FFT code was wrong.

[dpb]

OK, I did misunderstand what you meant there....

But you seem to be mistaking the overall time series envelope as necessarily indicating what individual frequency component's energy may be.

If you integrate over all frequency, you will return the total energy in the signal, but that's spread over all frequencies in the signal.

[dongwoo kim]

I wonder how the signal changes for each second rather than the FFT size of the whole signal. That is why the signal divided every 8 seconds was FFTed.

If you assume that the signal you divide in those eight seconds is the entire signal, the FFT value of the signal is still small.

[dpb]

The same effect is extant whatever subsection of the time series you look at...at least until you get so short a section that lower frequencies can't be sampled in that time.

The total energy is divided up amongst all frequencies so the integral over the full frequency band is the same. But, all frequencies are "mooshed together" at every instant of the time trace--there's no section that isn't made up of the summation of all frequencies no matter how short it is--even down to a single sample.

You seem to have some idea that the amplitude of a peak in the FFT should match the time signal amplitude--that's ONLY true for a pure, noise-free tone at exactly a frequency that matches a computed frequency bin identically. Even if the signal is pure, if the computed frequency bin isn't centered on that exact frequency the energy will be "smeared" across two or more bins and will still have to integrate to get the total energy.

[dongwoo kim]

Does that mean I have to change the time to sample low frequencies? How?

[dpb]

Back up and and slow down and think about what it is that's happening here...you've filtered the time domain signal between 130-150 Hz w/ a very strong filter so you've removed every bit of energy from the time domain signal outside that band before you did the FFT.

So, unless that original time series was also very much bandlimited, what is left of the original time series won't have anything at all close to the full timeseries magnitude to begin with, and then what energy is left is distributed amongst all the frequencies that do remain.

First, look at the time series of the filtered data and plot it, and compare it to the original.

[dongwoo kim]

When I plot, I wanted to compare the rising y-values of 130-150hz. So I filtered 130-150hz before doing FFT. Before that, I wanted to erase another high y-value. Because the Y-axis size of 130-150hz is too small.

I have compared according to your advice, but my problem is not solved.

[dpb]

What, precisely, do you think your problem really is?

What does the plot/magnitude of the filtered signal look like?

The y-axis amplitude is what it is owing to the two effects of having filtered out everything except 130-150 Hz signal to start with and then distributing that energy across the frequency bins making up that frequency range.

If you're looking for the total energy in that range, then perhaps it's the integral that you're looking for, not the individual peak value(s).

What, specifically, might be the answer is not possible to say because it isn't yet clear what the problem is that is trying to be solved.

[dongwoo kim]

My problem is that the signal with a large value has a small y-value added FFT. According to my assistant's hint, the peak value of the plot of that big signal should be large. I want to know the pressure signal changes over time (in seconds) through the FFT.

And the filtering was not f except 130-150 Hz, but FFT between 130-150 Hz?

[dpb]

Did you do as I suggested and plot the filtered time series? You didn't show it, if did.

I think your assistant may also be confused or he/she has a different idea of what to do than what you have done but I don't think there's any reason at all to expect a single peak to have a large value in the FFT for the reasons we've already been through (and, clearly, it doesn't).

"I want to know the pressure signal changes over time (in seconds) through the FFT. "

I don't know how to try to interpret the above...the FFT converts the time domain to the frequency domain; the assumption is the signal is stationary over the time frame of the transform. What the FFT (via the PSD you've calculated) shows you is the relative importance and magnitude of each frequency component (to the resolution of the number of points used) within the signal over the time period of the measurement. The time variable has been completely smeared in this representation just as the frequency domain is completely lost in the time domain.

"And the filtering was not f except 130-150 Hz, but FFT between 130-150 Hz?"

n=4; %% order of Butterworth filter
Wn=[130 150]; % 130 ~ 150Hz signal
Fn=Fs/2; % Nyquist
ftype='bandpass';
[b,a]=butter(n, Wn/Fn, ftype);%%% butter(): 버터워스 returns an order 2*n digital bandpass filter if
e=filtfilt(b,a,y);
N = length(e);
NFFT = 2^nextpow2(N); % Next power of 2 from length of y
Y = fft(e,NFFT);
...

I also don't know what the above was intended to say/ask, but you ran the original signal through a sharp (4 pole) bandpass filter between 130-150 Hz so you have effectively removed all the energy in the original signal outside that frequency range before doing the FFT/PSD. This is obvious if you look at the PSD plots--there's no energy outside that range.

Having done that, the magnitude of e is undoubtedly nothing like that of y which is, I believe, the time domain signal you're plotting and comparing magnitudes to. I'm guessing the magnitude of e ain't nothing close to that because you've filtered out most of it.

[dpb]

I've been trying to lead but don't seem to be getting anywhere...so finally decided to look for meself--

Fs = 1000; % Sampling frequency
T = 1/Fs;  % Sample time
L = 8000;  % Length of signal
t = (0:L-1)*T;                 % Time vector
n=4; %% order of Butterworth filter
Wn=[130 150]; % 130 ~ 150Hz signal
Fn=Fs/2; % Nyquist
ftype='bandpass';
[b,a]=butter(n, Wn/Fn, ftype);
y=TL26.x0_84708;
e=filtfilt(b,a,y);

reproduces your case far enough to look at both time series. So, what do we have???

>> [rms(y) rms(e)]
ans =
    ans =
    0.6610    0.0031
>>

which clearly shows that the magnitude of the filtered signal you're FFT'ing is NOT some large number but that your values are eminently reasonable.

figure
subplot(2,1,1)
ylabel('Original')
text(1,0.3,'RMS=0.661')
subplot(2,1,2)
plot(t,e)
ylabel('Filtered')
text(1,0.015,'RMS=0.0031)

results in

BTW, also not that

rms(detrend(y))
ans =
    0.1185

so a significant fraction of the magnitude of the signal is in the mean value, not the variability around the mean.

I don't have the advantage of knowing what the actual problem statement is here, but I'm guessing the FFT may not be the tool for whatever that task may be from what has been said. Perhaps simply the mean or rms of the signal over a series of time segments would be more nearly what are looking for--at least from the complaint about what this result is showing.

Again, without a precise problem description and statement of objective, it's virtually impossible to know the proper answer.

[dongwoo kim]

My experiment is to compare the pressure changes of a blocked pump.The unchanged graph measures the pressure of the pump when the tissue is not put in. The rest is to measure the pressure changes over time by putting in the tissue. I would like to see the pressure signal over time and know how much the amplitude and frequency of the frequency (BPF) for the tissues at runtime change. However, I am having a hard time because of the large amplitude of the non-waste vibrations.

The size of the FFT of the filtered file and the size of the unfiltered file are not very different, as I have compared the two signals by storing the amplitude every eight seconds.

Red point = constant pressure signal

Blue Point = Changing Pressure Signal

[dpb]

Well, maybe you're not measuring something that is sensitive to the change in conditions.

But, yes, in the specific energy band the two should be roughly the same altho the filter will have some reduction in amplitude owing to it.

I don't have any idea why you're using the bandpass filter centered around 140 Hz; the energy content in the signal itself is more nearly around 25 Hz. Below is overlaying the PSD for the same signal filtered and unfiltered. The center band of the filter is roughly same magnitude as the unfiltered; the rolloff on either side is apparent so there's no comparison outside that region.

But, the actual energy in the signal is way low in frequency comparatively; mostly what you're looking at in the 130-150 Hz band is probably just instrumentation noise.

As the linear scale shows, above 50 Hz there's nothing discernible at all in the unfiltered signal PSD; it only shows up in the filtered as a hump because you've totally cleaned out everything else to zero (to roughly machine precision).

Without a sample of the steady state and a transient with a known point at which the obstruction was introduced not much else could be done but I'm guessing a lowpass smoothing filter and just watching the rms value with time could possibly be the best option with these data.

Possibly redesigning the experiment may be the better route altho have no information on that but it doesn't seem as whatever is being measured/however it is being measured is particularly good for the purpose.

[dongwoo kim]

I compared the FFT graph with 'lab View' at the same time as I measured this pressure signal over time. The amplitude of where I want it in the lapview is very small. So I want to know how this part changes through Matlab's FFT. And frequenz is the number 140 by the calculation of the wings of the pump and RPM, and I know that the frequenz is changed by pressure, so I am curious about the value between 130-150 Hz.

The graph you showed above is time on the x-axis and amplitude on the y-axis?

[dpb]

"The graph you showed above is time on the x-axis and amplitude on the y-axis??"

No, it's the PSDs for the signal so the x-axis is frequency.

Whatever you calculate should be the beat frequency, you're either not measuring anything that reflects it or your sampling interval isn't correct in displaying the frequency axis or somesuch.

With the parameters you've given, there's essentially no energy at those frequencies in the measurement.

We have no access to anything but what you've posted, but if the attached data were sampled at 1 kHz, that's what your measurement shows. If it were actually being sampled at 10 kHz instead, that would look more nearly like what you're describing. But, we can't see your experimental setup from here.

if I instead use

Fs=10000;
f=Fs/2*linspace(0,1,NFFT/2+1);
plot(f,PSD1)
xlim([00 500])
ylim([0 0.06])

That's identically the same PSD, just making an assumption that the sampling frequency used before was in error. IF such were to be the case, then you could indeed see something in the frequency range you think you should.

But, the absolute peak is at 151 Hz, not 140 Hz altho there's energy in that area. Or, maybe the sampling frequency isn't precisely an even number?

What is the AC line frequency there? 50 Hz, maybe?

The upshot is "we can't tell" from the posted information alone.

[dongwoo kim]

How did you calculate the PSD1?

Can you show me the whole matlab code?

[dpb]

Same way you calculated Ya above except it is not using the filtered e series but the unfiltered y. It's the identical PSD as before, only the assumed sampling frequency was changed.

It is the identical data as in the lower of the two above plots "LinearPSD"(*) simply replotted by multiplying the f vector by 10 and then setting the range to 0-500. Or, you could get the same looking plot by just setting xlim([0 50]) on that plot--same data, just a change in assumed frequency scaling.

FFT() "knows nuthink!" about the sampling frequency; all it does is compute the amplitudes of each bin assuming a fixed sampling interval. The actual sampling rate and relationship to frequency is totally independent and based upon what the sampling rate was that collected the data. But, as you note, it is NOT part of the input to FFT() so can make no difference therein.

Can you verify for certain what the actual sampling rate was? Input a known fixed frequency sine wave into the same acquisition setup otherwise and see if you get the right output frequency.

What I did above was purely conjecture of a mistake that could possibly have been made -- no basis at all for doing so other than that is what would put the observed response energy at roughly the right location in frequency, again ASSUMING that that is also a correct statement of where that frequency should be. Again, we have no way to verify any of this.

(*) Actually, the blue line in both plots is the identical data; just the upper is on semilog scale and the lower is linear y scale as is the last above. But, the data is the same identical data in all three.

[dongwoo kim]

Actual sampling is unknown except for storing pressure signals 1000 times a second. So I set the sample rate at 1000.

I am investigating the pump BPF, and I know that this is changed by pressure signals. My assumption was that it would move between 130-150 Hz, and the actual value would have to be between 140 Hz So there can be a peak at 151hz.

[dpb]

"Actual sampling is unknown except for storing pressure signals 1000 times a second. So I set the sample rate at 1000. "

???? What does the above mean? You really are collecting data and don't know what the actual sampling rate is? That makes no sense at all.

I revert to the previous suggestion then as the only way I can see to proceed -- insert a known single-frequency sine wave into the measurement input in lieu of the pressure signal and figure out how to reproduce that with your data collection system.

Or, fix the system to record what you sample directly so it is known. Otherwise, looks to me as though you have insoluble conundrum.

[dongwoo kim]

I did not know that the pressure of the pump was stored 1000 times (1 khz) per second after the experiment. Then studying FFT required sampling at FFT, and I only did it at 1 khz.

Does this mean that the experimental sensor needs to be reset?

I have difficulty resetting the test sensor without knowing the details because my assistant has set up all the test sensors.

[dpb]

I have no idea what this means you must do...I have no way to know anything more about the setup nor data collection system.

All I know is that unless you can somehow confirm just what is actually being sampled at what rate and confirm you can get the right answer for a known input you have, essentially, nothing that is usable for any serious purpose.

Just scaling the input to match something you think should be is not satisfying at all.

Looks like you need to leave the office and go to the lab... :)

[dongwoo kim]

Thank you for all your advice so far :)

I'll go back to the lab based on your advice.

I wish you happiness in everything.

[dpb]

Good luck...glad to at least try to help! :)