Find the problem with my FFT code

1 次查看（过去 30 天）

显示更早的评论

dongwoo kim 2019-9-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/480461-find-the-problem-with-my-fft-code

评论： dpb 2019-9-17

在 MATLAB Online 中打开

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)|')

32 个评论
显示 30更早的评论隐藏 30更早的评论

dpb 2019-9-16

编辑：dpb 2019-9-16

在 MATLAB Online 中打开

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 2019-9-16

编辑：dpb 2019-9-16

在 MATLAB Online 中打开

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.

dpb 2019-9-17

编辑：dpb 2019-9-17

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.

dpb 2019-9-17

编辑：dpb 2019-9-17

在 MATLAB Online 中打开

"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.

dpb 2019-9-17

编辑：dpb 2019-9-17

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.