atan2演算法

⑸ 如何應用matlab進行fft分析

FFT是離散傅立葉變換的快速演算法，可以將一個信號變換
到頻域。有些信號在時域上是很難看出什麼特徵的，但是如
果變換到頻域之後，就很容易看出特徵了。這就是很多信號
分析採用FFT變換的原因。另外，FFT可以將一個信號的頻譜
提取出來，這在頻譜分析方面也是經常用的。
雖然很多人都知道FFT是什麼，可以用來做什麼，怎麼去
做，但是卻不知道FFT之後的結果是什意思、如何決定要使用
多少點來做FFT。

現在圈圈就根據實際經驗來說說FFT結果的具體物理意義。
一個模擬信號，經過ADC采樣之後，就變成了數字信號。采樣
定理告訴我們，采樣頻率要大於信號頻率的兩倍，這些我就
不在此羅嗦了。

采樣得到的數字信號，就可以做FFT變換了。N個采樣點，
經過FFT之後，就可以得到N個點的FFT結果。為了方便進行FFT
運算，通常N取2的整數次方。

假設采樣頻率為Fs，信號頻率F，采樣點數為N。那麼FFT
之後結果就是一個為N點的復數。每一個點就對應著一個頻率
點。這個點的模值，就是該頻率值下的幅度特性。具體跟原始
信號的幅度有什麼關系呢？假設原始信號的峰值為A，那麼FFT
的結果的每個點（除了第一個點直流分量之外）的模值就是A
的N/2倍。而第一個點就是直流分量，它的模值就是直流分量
的N倍。而每個點的相位呢，就是在該頻率下的信號的相位。
第一個點表示直流分量（即0Hz），而最後一個點N的再下一個
點（實際上這個點是不存在的，這里是假設的第N+1個點，也
可以看做是將第一個點分做兩半分，另一半移到最後）則表示
采樣頻率Fs，這中間被N-1個點平均分成N等份，每個點的頻率
依次增加。例如某點n所表示的頻率為：Fn=(n-1)*Fs/N。
由上面的公式可以看出，Fn所能分辨到頻率為為Fs/N，如果
采樣頻率Fs為1024Hz，采樣點數為1024點，則可以分辨到1Hz。
1024Hz的采樣率采樣1024點，剛好是1秒，也就是說，采樣1秒
時間的信號並做FFT，則結果可以分析到1Hz，如果采樣2秒時
間的信號並做FFT，則結果可以分析到0.5Hz。如果要提高頻率
分辨力，則必須增加采樣點數，也即采樣時間。頻率解析度和
采樣時間是倒數關系。
假設FFT之後某點n用復數a+bi表示，那麼這個復數的模就是
An=根號a*a+b*b，相位就是Pn=atan2(b,a)。根據以上的結果，
就可以計算出n點（n≠1，且n<=N/2）對應的信號的表達式為：
An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。
對於n=1點的信號，是直流分量，幅度即為A1/N。
由於FFT結果的對稱性，通常我們只使用前半部分的結果，
即小於采樣頻率一半的結果。

好了，說了半天，看著公式也暈，下面圈圈以一個實際的
信號來做說明。

假設我們有一個信號，它含有2V的直流分量，頻率為50Hz、
相位為-30度、幅度為3V的交流信號，以及一個頻率為75Hz、
相位為90度、幅度為1.5V的交流信號。用數學表達式就是如下：

S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

式中cos參數為弧度，所以-30度和90度要分別換算成弧度。
我們以256Hz的采樣率對這個信號進行采樣，總共采樣256點。
按照我們上面的分析，Fn=(n-1)*Fs/N，我們可以知道，每兩個
點之間的間距就是1Hz，第n個點的頻率就是n-1。我們的信號
有3個頻率：0Hz、50Hz、75Hz，應該分別在第1個點、第51個點、
第76個點上出現峰值，其它各點應該接近0。實際情況如何呢？
我們來看看FFT的結果的模值如圖所示。

圖1 FFT結果
從圖中我們可以看到，在第1點、第51點、和第76點附近有
比較大的值。我們分別將這三個點附近的數據拿上來細看：
1點： 512+0i
2點： -2.6195E-14 - 1.4162E-13i
3點： -2.8586E-14 - 1.1898E-13i

50點：-6.2076E-13 - 2.1713E-12i
51點：332.55 - 192i
52點：-1.6707E-12 - 1.5241E-12i

75點：-2.2199E-13 -1.0076E-12i
76點：3.4315E-12 + 192i
77點：-3.0263E-14 +7.5609E-13i

很明顯，1點、51點、76點的值都比較大，它附近的點值
都很小，可以認為是0，即在那些頻率點上的信號幅度為0。
接著，我們來計算各點的幅度值。分別計算這三個點的模值，
結果如下：
1點： 512
51點：384
76點：192
按照公式，可以計算出直流分量為：512/N=512/256=2；
50Hz信號的幅度為：384/(N/2)=384/(256/2)=3；75Hz信號的
幅度為192/(N/2)=192/(256/2)=1.5。可見，從頻譜分析出來
的幅度是正確的。
然後再來計算相位信息。直流信號沒有相位可言，不用管
它。先計算50Hz信號的相位，atan2(-192, 332.55)=-0.5236,
結果是弧度，換算為角度就是180*(-0.5236)/pi=-30.0001。再
計算75Hz信號的相位，atan2(192, 3.4315E-12)=1.5708弧度，
換算成角度就是180*1.5708/pi=90.0002。可見，相位也是對的。
根據FFT結果以及上面的分析計算，我們就可以寫出信號的表達
式了，它就是我們開始提供的信號。

總結：假設采樣頻率為Fs，采樣點數為N，做FFT之後，某
一點n（n從1開始）表示的頻率為：Fn=(n-1)*Fs/N；該點的模值
除以N/2就是對應該頻率下的信號的幅度（對於直流信號是除以
N）；該點的相位即是對應該頻率下的信號的相位。相位的計算
可用函數atan2(b,a)計算。atan2(b,a)是求坐標為(a,b)點的角
度值，范圍從-pi到pi。要精確到xHz，則需要采樣長度為1/x秒
的信號，並做FFT。要提高頻率解析度，就需要增加采樣點數，
這在一些實際的應用中是不現實的，需要在較短的時間內完成
分析。解決這個問題的方法有頻率細分法，比較簡單的方法是
采樣比較短時間的信號，然後在後面補充一定數量的0，使其長度
達到需要的點數，再做FFT，這在一定程度上能夠提高頻率分辨力。
具體的頻率細分法可參考相關文獻。

[附錄：本測試數據使用的matlab程序]
close all; %先關閉所有圖片
Adc=2; %直流分量幅度
A1=3; %頻率F1信號的幅度
A2=1.5; %頻率F2信號的幅度
F1=50; %信號1頻率(Hz)
F2=75; %信號2頻率(Hz)
Fs=256; %采樣頻率(Hz)
P1=-30; %信號1相位(度)
P2=90; %信號相位(度)
N=256; %采樣點數
t=[0:1/Fs:N/Fs]; %采樣時刻

%信號
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
%顯示原始信號
plot(S);
title('原始信號');

figure;
Y = fft(S,N); %做FFT變換
Ayy = (abs(Y)); %取模
plot(Ayy(1:N)); %顯示原始的FFT模值結果
title('FFT 模值');

figure;
Ayy=Ayy/(N/2); %換算成實際的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; %換算成實際的頻率值
plot(F(1:N/2),Ayy(1:N/2)); %顯示換算後的FFT模值結果
title('幅度-頻率曲線圖');

figure;
Pyy=[1:N/2];
for i="1:N/2"
Pyy(i)=phase(Y(i)); %計算相位
Pyy(i)=Pyy(i)*180/pi; %換算為角度
end;
plot(F(1:N/2),Pyy(1:N/2)); %顯示相點陣圖
title('相位-頻率曲線圖');

看完這個你就明白諧波分析了。

⑹ 如何使用opencv實現金字塔光流lk跟蹤演算法

#include <stdio.h>
#include <windows.h>
#include "cv.h"
#include "cxcore.h"
#include "highgui.h"
#include <opencv2\opencv.hpp>
using namespace cv;

static const double pi = 3.14159265358979323846;
inline static double square(int a)
{
return a * a;
}
/*該函數目的：給img分配內存空間，並設定format，如位深以及channel數*/
inline static void allocateOnDemand(IplImage **img, CvSize size, int depth, int channels)
{
if (*img != NULL) return;
*img = cvCreateImage(size, depth, channels);
if (*img == NULL)
{
fprintf(stderr, "Error: Couldn't allocate image. Out of memory?\n");
exit(-1);
}
}
/*主函數，原程序是讀取avi視頻文件，然後處理，我簡單改成從攝像頭直接讀取數據*/
int main(int argc, char *argv[])
{

//讀取攝像頭
VideoCapture cap(0);
//讀取視頻文件

//VideoCapture cap; cap.open("optical_flow_input.avi");
if (!cap.isOpened())
{
return -1;
}
Mat frame;

/*
bool stop = false;
while (!stop)
{
cap >> frame;
// cvtColor(frame, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);
// imshow("當前視頻", edges);
imshow("當前視頻", frame);
if (waitKey(30) >= 0)
stop = true;
}
*/

//CvCapture *input_video = cvCaptureFromFile( "optical_flow_input.avi" );
//cv::VideoCapture cap = *(cv::VideoCapture *) userdata;

//if (input_video == NULL)
// {
// fprintf(stderr, "Error: Can't open video device.\n");
// return -1;
// }

/*先讀取一幀，以便得到幀的屬性，如長、寬等*/
//cvQueryFrame(input_video);

/*讀取幀的屬性*/
CvSize frame_size;
frame_size.height = cap.get(CV_CAP_PROP_FRAME_HEIGHT);
frame_size.width = cap.get(CV_CAP_PROP_FRAME_WIDTH);

/*********************************************************/

/*用於把結果寫到文件中去,非必要
int frameW = frame_size.height; // 744 for firewire cameras
int frameH = frame_size.width; // 480 for firewire cameras
VideoWriter writer("VideoTest.avi", -1, 25.0, cvSize(frameW, frameH), true);

/*開始光流法*/
//VideoWriter writer("VideoTest.avi", CV_FOURCC('D', 'I', 'V', 'X'), 25.0, Size(640, 480), true);

while (true)
{
static IplImage *frame = NULL, *frame1 = NULL, *frame1_1C = NULL,
*frame2_1C = NULL, *eig_image = NULL, *temp_image = NULL,
*pyramid1 = NULL, *pyramid2 = NULL;

Mat framet;
/*獲取第一幀*/
// cap >> framet;
cap.read(framet);
Mat edges;
//黑白抽象濾鏡模式
// cvtColor(framet, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);

//轉換mat格式到lpiimage格式
frame = &IplImage(framet);
if (frame == NULL)
{
fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");
return -1;
}

/*由於opencv的光流函數處理的是8位的灰度圖，所以需要創建一個同樣格式的
IplImage的對象*/
allocateOnDemand(&frame1_1C, frame_size, IPL_DEPTH_8U, 1);

/* 把攝像頭圖像格式轉換成OpenCV慣常處理的圖像格式*/
cvConvertImage(frame, frame1_1C, 0);

/* 我們需要把具有全部顏色信息的原幀保存，以備最後在屏幕上顯示用*/
allocateOnDemand(&frame1, frame_size, IPL_DEPTH_8U, 3);
cvConvertImage(frame, frame1, 0);

/* 獲取第二幀 */
//cap >> framet;
cap.read(framet);
// cvtColor(framet, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);
frame = &IplImage(framet);
if (frame == NULL)
{
fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");
return -1;
}

/*原理同上*/
allocateOnDemand(&frame2_1C, frame_size, IPL_DEPTH_8U, 1);
cvConvertImage(frame, frame2_1C, 0);

/*********************************************************
開始shi-Tomasi演算法，該演算法主要用於feature selection,即一張圖中哪些是我
們感興趣需要跟蹤的點(interest point)
input:
* "frame1_1C" 輸入圖像.
* "eig_image" and "temp_image" 只是給該演算法提供可操作的內存區域.
* 第一個".01" 規定了特徵值的最小質量，因為該演算法要得到好的特徵點，哪就
需要一個選擇的閾值
* 第二個".01" 規定了像素之間最小的距離，用於減少運算復雜度，當然也一定
程度降低了跟蹤精度
* "NULL" 意味著處理整張圖片，當然你也可以指定一塊區域
output:
* "frame1_features" 將會包含fram1的特徵值
* "number_of_features" 將在該函數中自動填充上所找到特徵值的真實數目,
該值<= 400
**********************************************************/

/*開始准備該演算法需要的輸入*/

/* 給eig_image,temp_image分配空間*/
allocateOnDemand(&eig_image, frame_size, IPL_DEPTH_32F, 1);
allocateOnDemand(&temp_image, frame_size, IPL_DEPTH_32F, 1);

/* 定義存放frame1特徵值的數組，400隻是定義一個上限 */
CvPoint2D32f frame1_features[400];
int number_of_features = 400;

/*開始跑shi-tomasi函數*/
cvGoodFeaturesToTrack(frame1_1C, eig_image, temp_image,
frame1_features, &number_of_features, .01, .01, NULL);

/**********************************************************
開始金字塔Lucas Kanade光流法，該演算法主要用於feature tracking,即是算出
光流，並跟蹤目標。
input:
* "frame1_1C" 輸入圖像，即8位灰色的第一幀
* "frame2_1C" 第二幀，我們要在其上找出第一幀我們發現的特徵點在第二幀
的什麼位置
* "pyramid1" and "pyramid2" 是提供給該演算法可操作的內存區域，計算中間
數據
* "frame1_features" 由shi-tomasi演算法得到的第一幀的特徵點.
* "number_of_features" 第一幀特徵點的數目
* "optical_flow_termination_criteria" 該演算法中迭代終止的判別，這里是
epsilon<0.3，epsilon是兩幀中對應特徵窗口的光度之差的平方，這個以後的文
章會講
* "0" 這個我不知道啥意思，反正改成1就出不來光流了，就用作者原話解釋把
means disable enhancements. (For example, the second array isn't
pre-initialized with guesses.)
output:
* "frame2_features" 根據第一幀的特徵點，在第二幀上所找到的對應點
* "optical_flow_window" lucas-kanade光流演算法的運算窗口,具體lucas-kanade
會在下一篇詳述
* "5" 指示最大的金字塔層數，0表示只有一層，那就是沒用金字塔演算法
* "optical_flow_found_feature" 用於指示在第二幀中是否找到對應特徵值，
若找到，其值為非零
* "optical_flow_feature_error" 用於存放光流誤差
**********************************************************/

/*開始為pyramid lucas kanade光流演算法輸入做准備*/
CvPoint2D32f frame2_features[400];

/* 該數組相應位置的值為非零，如果frame1中的特徵值在frame2中找到 */
char optical_flow_found_feature[400];

/* 數組第i個元素表對應點光流誤差*/
float optical_flow_feature_error[400];

/*lucas-kanade光流法運算窗口,這里取3*3的窗口,可以嘗試下5*5,區別就是5*5
出現aperture problem的幾率較小,3*3運算量小，對於feature selection即shi-tomasi演算法來說足夠了*/
CvSize optical_flow_window = cvSize(5, 5);
// CvSize optical_flow_window = cvSize(5, 5);
/* 終止規則，當完成20次迭代或者當epsilon<=0.3，迭代終止，可以嘗試下別的值*/
CvTermCriteria optical_flow_termination_criteria= cvTermCriteria(CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3);

/*分配工作區域*/
allocateOnDemand(&pyramid1, frame_size, IPL_DEPTH_8U, 1);
allocateOnDemand(&pyramid2, frame_size, IPL_DEPTH_8U, 1);

/*開始跑該演算法*/
cvCalcOpticalFlowPyrLK(frame1_1C, frame2_1C, pyramid1, pyramid2,frame1_features, frame2_features, number_of_features,
optical_flow_window, 5, optical_flow_found_feature,optical_flow_feature_error, optical_flow_termination_criteria, 0);

/*畫光流場，畫圖是依據兩幀對應的特徵值，
這個特徵值就是圖像上我們感興趣的點，如邊緣上的點P(x,y)*/
for (int i = 0; i< number_of_features; i++)
{
/* 如果沒找到對應特徵點 */
if (optical_flow_found_feature[i] == 0)
continue;
int line_thickness;
line_thickness = 1;

/* CV_RGB(red, green, blue) is the red, green, and blue components
* of the color you want, each out of 255.
*/
CvScalar line_color;
line_color = CV_RGB(255, 0, 0);

/*畫箭頭,因為幀間的運動很小，所以需要縮放，不然看不見箭頭，縮放因子為3*/
CvPoint p, q;
p.x = (int)frame1_features[i].x;
p.y = (int)frame1_features[i].y;
q.x = (int)frame2_features[i].x;
q.y = (int)frame2_features[i].y;

double angle;
angle = atan2((double)p.y - q.y, (double)p.x - q.x);
double hypotenuse;
hypotenuse = sqrt(square(p.y - q.y) + square(p.x - q.x));

/*執行縮放*/
q.x = (int)(p.x - 5 * hypotenuse * cos(angle));
q.y = (int)(p.y - 5 * hypotenuse * sin(angle));

/*畫箭頭主線*/
/* "frame1"要在frame1上作畫.
* "p" 線的開始點.
* "q" 線的終止點.
* "CV_AA" 反鋸齒.
* "0" 沒有小數位.
*/
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);

/* 畫箭的頭部*/
p.x = (int)(q.x + 9 * cos(angle + pi / 4));
p.y = (int)(q.y + 9 * sin(angle + pi / 4));
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);
p.x = (int)(q.x + 9 * cos(angle - pi / 4));
p.y = (int)(q.y + 9 * sin(angle - pi / 4));
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);
}
/*顯示圖像*/

/*創建一個名為optical flow的窗口，大小自動改變*/
cvNamedWindow("Optical Flow", CV_WINDOW_NORMAL);
cvFlip(frame1, NULL, 2);
cvShowImage("Optical Flow", frame1);

/*延時，要不放不了*/
cvWaitKey(33);

/*寫入到文件中去*/

// cv::Mat m = cv::cvarrToMat(frame1);//轉換lpimgae到mat格式
// writer << m;//opencv3.0 version writer

}
cap.release();
cvWaitKey(33);
system("pause");
}

⑺ 在Flash AS3導彈跟蹤演算法里的，轉換為360度以內的正值不是很理解。

c2是轉換後的360度內的正值，則c2=(360-c1)%360-----此句有誤，轉換角度的方廳鉛法應該是：
c2=c1%360
................................應該是這樣吧？願共飢御同商討，扮肢好攜手進步...........張志晨..............

⑻ C語言迭代演算法初值問題

C語言沒有默認迭代陵棚初漏運值。
會自動初始化的有：
全局變數
靜態變數
在定義這兩種變數時，如果未賦初值尺搜則，則系統會將其值置為0.

你這里注釋部分有沒有執行的結果都一樣，可能是剛好你輸入的數據符合而已。
你再多試幾組數據看看。

⑼ 求這個碰撞演算法的解釋

肯定是一個公球和一個母球，呵呵

⑽ 求大神詳細講解c/c++/pascal凸包演算法

實這個演算法是在一年前得某場比賽中臨時抱佛腳學的，今天重新的來溫習了一遍
如何來理解凸包？一組平面上的點，求一個包含所有點的最小的凸多邊形，這就是凸包問題了。這可以形象地想成這樣：在地上放置一些不可移動的木樁，用一根繩子把他們盡量緊地圈起來，這就是凸包了，網路中的這張圖很生動+活潑+形象，所以你懂的

好說完這個我們首先要來了解下極角排序和左轉判定
極角排序：就是選取一個最左的點，按y最小，其次x最小來定義，接下來所有的點針對該點的射線，
按角度由小到大，若相同按距離由近到遠來排序
左轉判定：這個和叉積有關，對於向量p1(x1,y1),p2(x2,y2)如果x1*y2-x2*y1>0,則從p1到p2左轉

我學的是Graham演算法，那麼接下來來介紹下該演算法
（1）選取最下左的點P0
（2）計算出每個點相對於P0的角度和距離(利用這個來排序)排序
（3）設點數為n，將p[n-1]和p[0]入棧，判斷點集合是否為一條直線（初始k=2表示當前凸包的大小）
（4）i從1到n-1遍歷，對於p[k-1],p[k-2],p[i]若滿足左轉，將p[i]壓入棧
否則i--，k--
（5）k--,返回k表示凸包的點數

下面是我寫的模板

int Polygon::Graham(Polygon &con){//別用自己做對象
int t=0,i;
Point tmp;
//先y最小再x最小
for(i=1;i<n;i++)if(p[i]<p[t])t=i;
swap(p[t],p[0]);
for(i=0;i<n;i++){
tmp=p[i]-p[0];
p[i].dis=tmp.Len2();
p[i].angle=atan2(tmp.y,tmp.x);
}
sort(p,p+n,_cmp);
//for(int i=0;i<n;i++)p[i].out();
//cout<<"***"<<endl;
int k=0;
con.p[k++]=p[n-1];
con.p[k++]=p[0];
if(Sig(p[1].angle-p[n-1].angle)==0)con.p[k++]=p[n-1];//凸包為一線段
else{
for(i=1;i<n;i++){
//con[k-1].out();
//con[k-2].out();
//p[i].out();
if(Sig(Cross(con.p[k-1],con.p[k-2],p[i]))>0)con.p[k++]=p[i];
else {i--;k--;}
//cout<<"---"<<endl;
//for(int j=0;j<k;j++)con[j].out();
//system("pause");
}
}
return con.n=--k;
}
/*
9
1 4
3 6
5 7
2 2
3 3
5 4
8 3
9 6
7 1
*/

摘自http://blog.csdn.net/foreverlin1204/article/details/6221986