How to find corresponding frequencies for frequency magnitudes in FFT analysis

% I computed the fft for 10 values of xi and
% 10 values of l as shown below. I am trying to determine which values
% of xi and l oscillations occur. To achieve this,  tried to find the
% maximum freq. and use it as a threshold, so that if frequency < threshold
% freq. then find the xi and l values that correspond to the frequency <
% threshold. Please can anyone help with this?

clc
clear all

%%
mydata =load('sol_data.mat');


% x and y
x =mydata.x;
y =mydata.y;
nx = length(x);
ny = length(y);

% xi and l
xival=mydata.xival;
lval =mydata.lval;

% theta 
theta = mydata.theta;

X1 = mydata.X1;
Y1 = mydata.Y1;

% line in the HAN region
l = lval(end);

% x and y values at the centre of the channel
xm = (nx+1)/2;
ym = (ny+1)/2;

% Sampling frequency for planar and HAN regions
F = 1/(x(2)-x(1));     

% Length of x planar ad HAN regions
L = length(x);            

%% computing the fft
pfL = F/L*(0:L-1); 
fig2 = figure(2);
for ixi=1:length(xival)
    for il=1:length(lval)
        YplanarHAN(:, ixi, il)= abs(fft(theta(:,ixi, il)));
        hold on
        plot(pfL, YplanarHAN(:, ixi, il))
    end
end

% take the fftshift and exclude the zero values and remove the zero
% frequency
Fnotzero = [(-(L-1)/2):-1 1:(L-1)/2]*(F/L);
fig3 = figure(3);
for ixi=1:length(xival)
    for il=1:length(lval)

        % find the fft of theta
        YfftplanarHAN(:, ixi,il) = abs(fftshift(fft(theta(:,ixi, il))));

        % Remove the zero frequency 
        YfftnotzeroplanarHAN(:,ixi, il) = [YfftplanarHAN(1:((L-1)/2),ixi, il); YfftplanarHAN(((L+3)/2):L, ixi, il)];

        % plot
        hold on
        plot(Fnotzero, YfftnotzeroplanarHAN(:, ixi, il), 'LineWidth',1)
    end
end

%% Calculate the peak frequency

% length of fftshift of the FFT of the theta
LFnotzero = 1:numel(Fnotzero);           
% Calculate the peak frequency
for ixi=1:length(xival)
    for il=1:length(lval)
        [peak_mag(:, ixi, il), locs(:, ixi, il)] =  findpeaks(abs(YfftnotzeroplanarHAN(LFnotzero)), 'MinPeakProminence',0.25);
    end
end 

% Frequency corresponding to the peak magnitudes
Freq_peak = Fnotzero(locs);

% Maximum frequency for each xi 
Max_freq = max(Freq_peak(:));

% Threshold frequency
threshold_max = max(Freq_peak(:));

% Test
test1 = abs(Fnotzero) < threshold_max;
[xivals1, Lvals1] = find(test1);
xi_1 = xival(xivals1);
Lv_1 = lval(Lvals1);

Index exceeds the number of array elements. Index must not exceed 10.

Answer 1

[Star Strider]

I am not certain what you want from these calculations.

I used my own fft call (since I found it difficult to understand your code), and then used findpeaks to find the peaks and their frequencies and then returned them in cell arrays ‘Peaks’ and ‘Freqs’ with the second subscript (‘k2’) denoting the ‘lval’ subscript value (not the actual ‘lval’ value) and the first subscript (‘k1’) denoting the ‘page’ (third dimension) of the ‘FTtheta’ (Fourier Transform of ‘theta’) array. So the subscripts are essentially {page subscript, ‘lval’ subscript) for each element of each cell array.

% 10 values of l as shown below. I am trying to determine which values
% of xi and l oscillations occur. To achieve this,  tried to find the
% maximum freq. and use it as a threshold, so that if frequency < threshold
% freq. then find the xi and l values that correspond to the frequency <
% threshold. Please can anyone help with this?

% clc
% clear all
% 
% %%
mydata =load('sol_data.mat');


% x and y
x =mydata.x;
y =mydata.y;
nx = length(x);
ny = length(y);

% xi and l
xival=mydata.xival;
lval =mydata.lval;

% theta 
theta = mydata.theta;

X1 = mydata.X1;
Y1 = mydata.Y1;

% line in the HAN region
l = lval(end);

% x and y values at the centre of the channel
xm = (nx+1)/2;
ym = (ny+1)/2;

% Sampling frequency for planar and HAN regions
F = 1/(x(2)-x(1));     

% Length of x planar ad HAN regions
L = length(x);            

FTtheta  = fft(theta .* hann(size(theta,1)))/size(theta,1);
FTtheta = fftshift(FTtheta,1);
Fv = x;
Iv = 1:numel(Fv);

figure
hold on
xx = repmat(Fv,numel(lval),1).';
yy = repmat(lval,numel(Fv),1);
for k = 1:size(FTtheta,3)
    plot3(xx, yy, abs(FTtheta(Iv,:,k)))
end
hold off
grid
xlabel('Frequency')
ylabel('lval')
view(-27,30)

figure
for k1 = 1:size(FTtheta,3)
    figure
    plot3(xx, yy, abs(FTtheta(Iv,:,k1)))
    grid on
    xlabel('Frequency')
    ylabel('lval')
    title("Page "+k1)
    for k2 = 1:size(FTtheta,2)
        [pks, frqs] = findpeaks(abs(FTtheta(Iv,k2,k1)), Fv, 'MinPeakProminence',0.01);
        Peaks{k1,k2} = pks;
        Freqs{k1,k2} = frqs;
    end
    xlim([-1 1]*0.2E-5)
end

% Illustrating The Cell Array — 

First_Peaks_Cell = Peaks{1,1}

First_Peaks_Cell = 2×1

0.1574 0.1574

First_Freqs_Cell = Freqs{1,1}

First_Freqs_Cell = 1×2

1.0e-06 * -0.7333 0.7333

First_lval_Value = lval(1)

First_lval_Value = 2.0000e-08

% ----------
Second_Peaks_Cell = Peaks{1,2}

Second_Peaks_Cell = 2×1

0.0724 0.0724

Second_Freqs_Cell = Freqs{1,2}

Second_Freqs_Cell = 1×2

1.0e-06 * -0.7333 0.7333

Second_lval_Value = lval(2)

Second_lval_Value = 7.5852e-07

% ----------
Middle_Peaks_Cell = Peaks{5,5}

Middle_Peaks_Cell = 0.0933

Middle_Freqs_Cell = Freqs{5,5}

Middle_Freqs_Cell = 0

Middle_lval_Value = lval(5)

Middle_lval_Value = 2.9741e-06

% ----------
Last_Peaks_Cell = Peaks{10,10}

Last_Peaks_Cell = 0.2414

Last_Freqs_Cell = Freqs{10,10}

Last_Freqs_Cell = 0

Last_lval_Value = lval(10)

Last_lval_Value = 6.6667e-06

return                                      % Stop Here

%% computing the fft
pfL = F/L*(0:L-1); 
fig2 = figure(2);
for ixi=1:length(xival)
    for il=1:length(lval)
        YplanarHAN(:, ixi, il)= abs(fft(theta(:,ixi, il)));
        hold on
        plot(pfL, YplanarHAN(:, ixi, il))
    end
end

% take the fftshift and exclude the zero values and remove the zero
% frequency
Fnotzero = [(-(L-1)/2):-1 1:(L-1)/2]*(F/L);
fig3 = figure(3);
for ixi=1:length(xival)
    for il=1:length(lval)

        % find the fft of theta
        YfftplanarHAN(:, ixi,il) = abs(fftshift(fft(theta(:,ixi, il))));

        % Remove the zero frequency 
        YfftnotzeroplanarHAN(:,ixi, il) = [YfftplanarHAN(1:((L-1)/2),ixi, il); YfftplanarHAN(((L+3)/2):L, ixi, il)];

        % plot
        hold on
        plot(Fnotzero, YfftnotzeroplanarHAN(:, ixi, il), 'LineWidth',1)
    end
end

%% Calculate the peak frequency

% length of fftshift of the FFT of the theta
LFnotzero = 1:numel(Fnotzero);           
% Calculate the peak frequency
for ixi=1:length(xival)
    for il=1:length(lval)
        [peak_mag(:, ixi, il), locs(:, ixi, il)] =  findpeaks(abs(YfftnotzeroplanarHAN(LFnotzero)), 'MinPeakProminence',0.25);
    end
end 

% Frequency corresponding to the peak magnitudes
Freq_peak = Fnotzero(locs);

% Maximum frequency for each xi 
Max_freq = max(Freq_peak(:));

% Threshold frequency
threshold_max = max(Freq_peak(:));

% Test
test1 = abs(Fnotzero) < threshold_max;
[xivals1, Lvals1] = find(test1);
xi_1 = xival(xivals1);
Lv_1 = lval(Lvals1);

.

[University]

Thank you so much. I will check your code in detail to understand what you have done. Essentially, I want determine where there oscillations in the system. The idea is that after finding the peak frequency, I want to set a threshold frequency such that if the frequency < threshold frequency, it means no oscillation and I want to return the xivals and lvals for such frequencies.

[Star Strider]

My pleasure!

I do not entirely understand what you want to do. My code returns the peak values and their frequencies (I believe those are the ‘xvals’) and the index into the ‘lval’ vector, so it should provide the essential information for you to process the data further to determine the threshold frequencies.