Getting amplitude from a noisy signal

I'm working on analyzing some data from a grain boundary in an atomic simulation; the included file has the relevant variables. I've got a lot of somewhat similar data to go through and am trying to determine a reliable way to get the amplitude of my data. There are a few complications that make the different tools I'm using difficult:

1. The x-axis location sometimes shifts forward and backwards, but usually by negligible amounts. My included code sorts it to be ascending simply because many solutions require it, but an implementation that doesn't need to sort would reduce error.
2. The data has noise that make a lot of "false peaks" and "false valleys" that play havoc on using a simple findpeak function. Unfortunately, the sample rate for my data isn't high enough to make this easy. Perhaps leveraging the periodic behavior of the signal could help?
3. My attempts to smooth out the noise with a filter drastically reduce the amplitude, when that is the variable I am trying to study. I've been using a 2nd order savitzky-golay filter via the sgolayfilt() command, but I do not know if that is the best fit for this work.

I have very little experience with signal processing, and only somewhat more with MatLab. I think a fourier decomposition and then reconstructing the signal would be my next best option, but after plotting the power spectrum and deciding that somewhere around 0.36 angstroms would be a good cutoff I'm not sure how to proceed. I feel like this would get an accurate periodicity, and then perhaps I could use a combination of movmean with movmax and movmin to get the original amplitude. It isn't yet clear if the periodicity varies slightly.

The following code plots the power spectrum and frequency, then the actual data (sorted along x-axis).

%% Load data
Temp = '10';
fileName = ['powerDump', Temp, '.mat'];
load(fileName);


%% Sort data by position

[~,sortID] = sort(newGBLocation);

%% compute power spectrum

% pspectrum expects increasing second argument, but it moves erratically in
% the negative Z-direction and needs sorting
[p,f] = pspectrum(newPEAverage(sortID),newGBLocation(sortID)-min(newGBLocation));

% plot energy vs. position data [Angstroms]
figure
plot(newGBLocation(sortID)-min(newGBLocation),newPEAverage(sortID),'r-')
xlabel('Position $[\AA]$','interpreter','latex')
ylabel('Potential Energy [eV]','interpreter','latex')

% plot power vs. wavelength [Angstroms]
figure; 
plot(1./f,p);
xlabel('$\lambda [\AA]$','interpreter','latex')
ylabel('Power','interpreter','latex')
set(gca,'yscale','log','xscale','log')

% plot power spectrum (power vs. frequency [1/Angstroms])
figure; 
plot(f,p); 
xlabel('$f [\AA^{-1}]$','interpreter','latex')
ylabel('Power','interpreter','latex')
set(gca,'yscale','log','xscale','log')

[Pat Gipper]

Is the area I circled what you consider "signal"?

If so, then you could feed the re-ordered samples into a bandpass filter centered around 0.36 Angstroms. I used a 4th order filter with a cut-in wavelength of 1 A and stop wavelength of 0.0125 A and that seemed to work. I did remove the mean value of -4.41 eV first before feeding that data into the filter.

[Eric]

This seems really promising! What function/app did you use for this? bandpass() doesn't seem to have a way to specify the order of filtering, and I can't find much about a cut-in wavelength.

[Alan]

Hi Eric,

You can use butter() function to generate filter parameters for a 4th order bandpass filter, and then apply the filter by using the sosfilt() function:

[z,p,k] = butter(4, [f1 f2]/(fs/2), 'bandpass'); % Where f1 and f2 are the frequency limits of the pass band

[sos, g] = zp2sos(z,p,k);

new_newPEAverage = sosfilt(sos, newPEAverage);

I’m assuming that the irregularities in newGBLocation that you mention are the tiny bumps in the plot below:

If the points in newGBLocation do not directly map to points in newPEAverage, and if we can assume that the points in newGBLocation are supposed to represent a straight line, polyfit() can be used to fit a line:

len = size(newGBLocation, 2);

p = polyfit(1:len, newGBLocation, 1);

new_newGBLocation = polyval(p, 1:len);

The following is the plot of new_newGBLocation:

Otherwise, if the data points in newGBLocations can be directly mapped to newPEAverage, your approach of sorting the data points should do the job.

Here are documentation pages for the functions I’ve used:

1. butter() : https://www.mathworks.com/help/signal/ref/butter.html
2. zp2sos() : https://www.mathworks.com/help/signal/ref/zp2sos.html
3. sosfilt() : https://www.mathworks.com/help/signal/ref/sosfilt.html
4. polyfit() : https://www.mathworks.com/help/matlab/ref/polyfit.html
5. polyval() : Polynomial evaluation - MATLAB polyval - MathWorks India

There are various filters other than the Butterworth filter mentioned in this page: https://www.mathworks.com/help/signal/ug/iir-filter-design.html

I hope this helps.