How to find corresponding frequencies for frequency magnitudes in FFT analysis

4 views (last 30 days)
% 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.

Accepted Answer

Star Strider
Star Strider on 27 Feb 2024
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);
.
  2 Comments
University
University on 27 Feb 2024
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
Star Strider on 27 Feb 2024
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Time-Frequency Analysis in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!