Find Periodicity of a plot to do Fourier Transformation

1 回表示 (過去 30 日間)
Alex Perrakis
Alex Perrakis 2021 年 11 月 22 日
コメント済み: Star Strider 2021 年 11 月 22 日
Hello Guys, i have following code, and which produces some plots, i would like to find the periodicity of the plot to do a Fourier. But it does not work, maybe one of you can tell me why?
pic01=imread("1.tiff");
pic050=im2gray( imread("50.tiff") );
newpicture=imsubtract(pic01,pic050);
%imshow(newpicture)
% level = graythresh(newpicture);
level = 0.1405;
newbinpic = imbinarize(newpicture,level);
%imshowpair(newpicture,newbinpic,'montage')
newbinpic2=bwpropfilt(newbinpic,'perimeter',1);
% imshowpair(newbinpic2,newpicture,'montage');
boundaries = bwboundaries(newbinpic2,'noholes');
% Plot raw (x,y) values
x = boundaries{1}(:,2);
y = boundaries{1}(:,1);
figure()
subplot(2,2,1)
plot(x,y, 'r-')
title('raw values')
axis equal tight
hold on
% Assuming (x,y) is a noisy circle, find center
cnt = [min(x)+range(x)/2, min(y)+range(y)/2];
xline(cnt(1))
yline(cnt(2))
% Shift to (0,0) and compute polar coordinates
[theta, rho] = cart2pol(x-cnt(1), y-cnt(2));
% Plot in polar coord
subplot(2,2,2)
polarplot(theta, rho)
title('Polar coord.')
% Define a reasonable threshold to isolate the spike in RHO values
% 90% of the median looks OK
subplot(2,2,3)
plot(rho)
threshold = median(rho)*.9;
yline(threshold, '-k', 'Threshold')
ylabel('rho'); xlabel('index')
title('rho values')
% Find start and end of peak
peakIdx = find(rho < threshold);
peakIdx(2:end-1) = [];
xline(peakIdx(1)) % plot start of peak
xline(peakIdx(2)) % plot end of peak
% Remove the peak
theta(peakIdx(1):peakIdx(2)) = []; % you can do the same with x
rho(peakIdx(1):peakIdx(2)) = []; % you can do the same with y
% Plot the circle without the interior object
subplot(2,2,4)
plot(theta,rho)
t=[1:1:1730]';
title('Inner removed')
[pks,locs] = findpeaks(rho(:,1)); % Peaks & Locations
mean_period = mean(diff(t(locs))); % Period
text(min(thetalim)+0.1*diff(thetalim), min(thetalim)+0.95*diff(thetalim), sprintf('Mean Period = %.2f time units', mean_period))
I attached the Pictures, but because it does not accept .tiff , i uploaded in png, just change it in the code.
  8 件のコメント
Alex Perrakis
Alex Perrakis 2021 年 11 月 22 日
And how am i able to compute this frequency? As you can the only data i have is those two images.
Star Strider
Star Strider 2021 年 11 月 22 日
For spatial data, the sampling interval would be the pixel interval in each dimension. Add spatial calibration (e.g. parsecs/pixel) if available and necessary. The rest follows directly.
.

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangePolar Plots についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by