0

I have collected some data of pressure distribution of an acoustic field along an axis, i.e. Pressure_vs_distance, at a specific source frequency (Let's say 0.5 MHz). Now, I want to obtain the frequency components of this data (which will be n*0.5 MHz) but every code I try, I can't extract those frequencies and the FFT diagram is like a vertical line at zero. The data and figure are provided below. Could you please help me to see what's wrong?

P.S: My own guess is the number of axial data. But, even when I increased it, the result didn't change (shape of the fft diagram).

Pressure (MPa) _vs_axial distance (cm)

data (Size~ 2 Mb): https://ufile.io/wjhlo

L=length(y);
dx=9.43479300064157e-05;
fs=1/dx;
out=fft(y,L)/L;
figure
plot(fs/2*linspace(0,1,(length(out)/2)+1),abs(out(1:(length(out)/2)+1))) 
title('One sided Spectrum')
xlabel('Normalized frequency')
ylabel('Magnitude')
Farhad
  • 3
  • 1
  • 2
    This is because the 0th frequency is much, much stronger than any other frequency. The 0th frequency is the DC component. Subtracting the DC component (`y-mean(y)`) or deleting the 0th frequency component (`out(1)=0`) will allow you to see your plot the way you expect it. – Cris Luengo May 17 '18 at 14:51
  • Thank you for your response. I tried that, but it didn't work. It's not exactly and only at zero. If you could test the data yourself, that'd be really helpful. – Farhad May 17 '18 at 15:07
  • Then change the y limits: `set(gca,'ylim',[0,10])`, replace `10` until you see what you want to see. – Cris Luengo May 17 '18 at 15:09
  • Or use a log-lin plot: `set(gca,'yscale','log')`. – Cris Luengo May 17 '18 at 15:10
  • The problem is I can't see the peaks at expected frequencies, 0.5*n MHz (n=1,2,...) (x-axis). Actually, it's a flat line at high frequencies (y=0 @ f>20 Hz) – Farhad May 17 '18 at 15:20

2 Answers2

2

From this question:

Confusion in figuring out the relation between actual frequency values and FFT plot indexes in MATLAB

[Ycomp, fHz] = getFFT(y(1:1000),1/(x(2)-x(1)))

enter image description here

Your sampling frequency seems to be 10.6KHz. So according to the Nyquist Shannon sampling theorem you will not detect frequencies above 5.3 KHz. Not sure where your 0.5 MHz is coming from.

Gelliant
  • 1,835
  • 1
  • 11
  • 23
1

First of all prepare your data: try pass your data through the pre-emphasis filter:

y1 = filter([1 -1], 1, y);

Next, observe your data and remove obvious data errors:

y1(1:3) = 0;

The Fourier analysis suites only for quasi-periodic signals. So, you should select an interval of quasi-stationarity and perform analysis on such (possible overlapping) intervals. Let's say the signal is quasi-stationary on 5000 points. And we want to perform analysis with 100 points shift:

sampleRate = 1/(x(2)-x(1));
frameSize = 5000;
frameShift = 100;
[signalSpectrogram, freqObs, timeObs] = spectrogram(y1, frameSize, frameSize-frameShift, pow2(nextpow2(frameSize)), sampleRate);

Let's display all the results:

figure('NumberTitle','off', 'Name','The Spectral Analysis', 'Units','normalized', 'Position',[0 0 1 1]);

subplot(3,1,1);
plot(x,y);
grid('on');
xLim = xlim();
title('Original Signal');

subplot(3,1,2);
plot(x,y1);
grid('on');
title('Prepared Signal');

subplot(3,1,3);
imagesc(timeObs, freqObs, 10*log10(signalSpectrogram.*conj(signalSpectrogram)));
axis('xy');
axis([xLim 0 200]);
title('Prepared Signal Spectrogram');

At the end you must get something similar to: spectrogram

Andrei Davydov
  • 315
  • 1
  • 7