原理
短时傅里叶变换(Short Time Fourier Transform, STFT) 是一个用于语音信号处理的通用工具.它定义了一个非常有用的时间和频率分布类, 其指定了任意信号随时间和频率变化的复数幅度. 实际上,计算短时傅里叶变换的过程是把一个较长的时间信号分成相同长度的更短的段, 在每个更短的段上计算傅里叶变换, 即傅里叶频谱.
短时傅里叶变换通常的数学定义如下:
其中,
DTFT (Decrete Time Fourier Transform) 为离散时间傅里叶变换.  其数学公式, 如下所示:
其中,  x(n) 为在采样数 n 处的信号幅度. ω~ 的定义如下:
实现时, 短时傅里叶变换被计算为一系列加窗数据帧的快速傅里叶变换 (Fast Fourier Transform, FFT)，其中窗口随时间 “滑动” (slide) 或“跳跃” (hop) 。
Python 实现
在程序中, frame_size 为将信号分为较短的帧的大小, 在语音处理中, 通常帧大小在 20ms 到 40ms 之间. 这里设置为 25ms, 即 frame_size = 0.025;
frame_stride 为相邻帧的滑动尺寸或跳跃尺寸, 通常帧的滑动尺寸在 10ms 到 20ms 之间, 这里设置为 10ms, 即 frame_stride = 0.01. 此时, 相邻帧的交叠大小为 15ms;
窗函数采用汉明窗函数 (Hamming Function) ;
在每一帧, 进行 512 点快速傅里叶变换, 即 NFFT = 512. 具体程序如下:
# -*- coding: utf8 -*-
import numpy as np
def calc_stft(signal, sample_rate=16000, frame_size=0.025, frame_stride=0.01, winfunc=np.hamming, NFFT=512):
# Calculate the number of frames from the signal
frame_length = frame_size * sample_rate
frame_step = frame_stride * sample_rate
signal_length = len(signal)
frame_length = int(round(frame_length))
frame_step = int(round(frame_step))
num_frames = 1 + int(np.ceil(float(np.abs(signal_length - frame_length)) / frame_step))
# zero padding
pad_signal_length = num_frames * frame_step + frame_length
z = np.zeros((pad_signal_length - signal_length))
# Pad signal to make sure that all frames have equal number of samples
# without truncating any samples from the original signal
pad_signal = np.append(signal, z)
# Slice the signal into frames from indices
indices = np.tile(np.arange(0, frame_length), (num_frames, 1)) + \
np.tile(np.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
frames = pad_signal[indices.astype(np.int32, copy=False)]
# Get windowed frames
frames *= winfunc(frame_length)
# Compute the one-dimensional n-point discrete Fourier Transform(DFT) of
# a real-valued array by means of an efficient algorithm called Fast Fourier Transform (FFT)
mag_frames = np.absolute(np.fft.rfft(frames, NFFT))
# Compute power spectrum
pow_frames = (1.0 / NFFT) * ((mag_frames) ** 2)
return pow_frames
if __name__ == '__main__':
import scipy.io.wavfile
import matplotlib.pyplot as plt
# Read wav file
# "OSR_us_000_0010_8k.wav" is downloaded from http://www.voiptroubleshooter.com/open_speech/american.html
sample_rate, signal = scipy.io.wavfile.read("OSR_us_000_0010_8k.wav")
# Get speech data in the first 2 seconds
signal = signal[0:int(2. * sample_rate)]
# Calculate the short time fourier transform
pow_spec = calc_stft(signal, sample_rate)
plt.imshow(pow_spec)
plt.tight_layout()
plt.show()
参考资料
1. DISCRETE TIME FOURIER TRANSFORM (DTFT). https://www.dsprelated.com/freebooks/mdft/Discrete_Time_Fourier_Transform.html
2. THE SHORT-TIME FOURIER TRANSFORM. https://www.dsprelated.com/freebooks/sasp/Short_Time_Fourier_Transform.html
3. Short-time Fourier transform. https://en.wikipedia.org/wiki/Short-time_Fourier_transform
4. Speech Processing for Machine Learning: Filter banks, Mel-Frequency Cepstral Coefficients (MFCCs) and What's In-Between. https://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html
Python scipy 计算短时傅里叶变换（Short-time Fourier transforms)
计算短时傅里叶变换(STFT) scipy.signal.stft(x,fs = 1.0,window ='hann',nperseg = 256,noverlap = None,nfft = Non ...
matlab 时频分析（短时傅里叶变换、STFT）
短时傅里叶变换,short-time fourier transformation,有时也叫加窗傅里叶变换,时间窗口使得信号只在某一小区间内有效,这就避免了传统的傅里叶变换在时频局部表达能力上的不足, ...
傅里叶变换 - Fourier Transform
傅里叶级数 傅里叶在他的专著<热的解析理论>中提出,任何一个周期函数都可以表示为若干个正弦函数的和,即: \[f(t)=a_0+\sum_{n=1}^{\infty}(a_ncos(n\o ...
「学习笔记」Fast Fourier Transform
前言 快速傅里叶变换(\(\text{Fast Fourier Transform,FFT}\) )是一种能在\(O(n \log n)\)的时间内完成多项式乘法的算法,在\(OI\)中的应用很多,是 ...
浅谈范德蒙德&lpar;Vandermonde&rpar;方阵的逆矩阵的求法以及快速傅里叶变换&lpar;FFT&rpar;中IDFT的原理
浅谈范德蒙德(Vandermonde)方阵的逆矩阵与拉格朗日(Lagrange)插值的关系以及快速傅里叶变换(FFT)中IDFT的原理 标签: 行列式 矩阵 线性代数 FFT 拉格朗日插值 只要稍微看 ...
短网址&lpar;short URL&rpar;系统的原理及其实现
短网址(short URL)系统的原理及其实现 https://hufangyun.com/2017/short-url/?hmsr=toutiao.io&utm_medium=toutiao ...
paip&period;编程语言方法重载实现的原理及python&comma;php&comma;js中实现方法重载
paip.编程语言方法重载实现的原理及python,php,js中实现方法重载 有些语言,在方法的重载上,形式上不支持函数重载,但可以通过模拟实现.. 主要原理:根据参数个数进行重载,或者使用默认值  ...
MapReduce 原理与 Python 实践
MapReduce 原理与 Python 实践 1. MapReduce 原理 以下是个人在MongoDB和Redis实际应用中总结的Map-Reduce的理解 Hadoop 的 MapReduce  ...
随机推荐
&lbrack;C&num;&rsqb;&period;NET中几种Timer的使用
这篇博客将梳理一下.NET中4个Timer类,及其用法. 1. System.Threading.Timer public Timer(TimerCallback callback, object s ...
&lbrack;BZOJ1131&rsqb;&lbrack;POI2008&rsqb; Sta 树的深度
Description 给出一个N个点的树,找出一个点来,以这个点为根的树时,所有点的深度之和最大 Input 给出一个数字N,代表有N个点.N<=1000000 下面N-1条边. Output ...
C&num;&lowbar;简单实用的翻页
简单实用的生成翻页HTML辅助类 C# using System.Text; namespace ClassLibrary { ///  /// /// ...
Equivalent Strings
Equivalent Strings 题目链接:http://acm.hust.edu.cn/vjudge/contest/view.action?cid=84562#problem/E 题意: 给出 ...
IIS-Server is too busy &lowbar;解决方法
httpRuntime Server Too Busy 修改方法:修改服务器.net配置“machine.config＂文件,该文件位于Windows系统目录下,如“C:\WINDOWS \Micro ...
（译）详解javascript立即执行函数表达式（IIFE）
写在前面 这是一篇译文,原文:Immediately-Invoked Function Expression (IIFE) 原文是一篇很经典的讲解IIFE的文章,很适合收藏.本文虽然是译文,但是直译的 ...
ubuntu启动慢
http://blog.sina.com.cn/s/blog_4cc9ffbc0100rxhh.html 参考 笔记本装的是ubuntu12.04,最近发现开机启动特别慢,至少3分钟,总担心系统会启动 ...
POJ 1325、ZOJ 1364、HDU 1150 Machine Schedule - from lanshui&lowbar;Yang
Problem Description As we all know, machine scheduling is a very classical problem in computer scien ...
一个简单的redis调用类
能只能判断函数的调用规则,容错规则, 例如set函数 set($key, $value, $time = false) 根据time的真假来判断是否使用set,或者是setex函数 get函数 get ...
oracle修改字符集后数据库不能启动
最近在做修改字符集的实验,悲剧的是修改后重启,数据库启动不了. SQL> alter system set nls_language='AMERICAN'   scope = spfile; S ...

lemon lelieven：MRI原理-信号​zhuanlan.zhihu.com
之前已经得到了MRI的信号公式：
不过这是一个连续的信号，经过A/D转换、采样之类的操作，变成离散信号，才能使用计算机处理。
ADC的过程是在频率编码梯度打开时进行采集的，每个采集点的值是：
变换一下就是：
可见采样就直接得到了k空间。
每进行一次相位编码，采样一次，得到一行k空间数据，重复这个过程就可以得到足够的数据来重建图像。
得到K空间之后逆变换就可以得到图像。
二维的DFT变换对是：
从采样开始，一直到最后的图像重建，基本都是信号处理和数字图像处理的内容，这里不再细说，知乎里有各位大佬的全方位科普。比如，
阿姆斯特朗：形象理解二维傅里叶变换​zhuanlan.zhihu.com
再废话几句。
实际中傅里叶变换的信号有以下几种：
（1） 时域：有限非周期连续信号，频域：连续信号
（2） 时域：周期连续信号，频域：离散、无限信号
（3） 时域：有限、离散信号，频域：连续、周期信号
（4） 时域：周期、离散信号，频域：离散、周期信号
第四种情况下，时域、频域只取其中一个周期处理，就是实际处理中的离散傅里叶变换（DFT）。
实际中，K空间存储的都是复数，二维K空间相当于一个二维复空间，相当于一个4维实空间，实际中是无法画出来的。
上面的图，大多书告诉你，左边就是K空间。实际左边只是K空间的幅度图，再加上相位图才是完整的K空间。
下面来几张动图直接感受一下二维傅里叶变换的情况。
K空间只有一个点时，对应的图像就是一个明暗变化的条纹，是一个非常单纯的二维波形。中心代表频率为0，往周边频率增高，右侧条纹周期降低，频率增高。并且条纹的方向与左侧点的方向垂直
上面三幅说明了K空间变化对图像的影响。当数据从中心往两边增加时，可以发现先出现图像的亮度，随后才会有完整的组织对比度。，说明低频信息影响图像的轮廓，高频信息影响图像的对比度。从傅里叶变换的频率分解角度来看，也是一致的。
最后说一个国内医学影像界对K空间和傅里叶变换流传很广的一个错误。
最早可能来自下面这本书。
这是书里第47页的内容，这一大段，明显全错，作者明显不明白傅里叶变换到底时怎么一回事。二维傅里叶变换是两个频率，相位编码只是相当于其中的一个频率。相位编码的次数和频率编码的次数时相互独立的，没有关系。如果继续看这本书后面对K空间的描述，会发现作者把傅里叶变换的相位、成像过程中的相位编码明显当作相位来看，把频率编码当作频率来看。如果有人看过这本书的配套视频，作者似乎也确实时这么认为的，可惜时错的啊。
本书2008年发行，又是一众国内大佬联合编写，发行量巨大，错误流传也很广。时不时看到有人问为什么傅里叶变换只能区分180°的相位差？当然时作者错了啊，还能为什么。
最后个人经验来看，不经过手算傅里叶积分，没有亲自算过图像的傅里叶变换，很难说真正懂了K空间和傅里叶变换。

原理
短时傅里叶变换(Short Time Fourier Transform, STFT) 是一个用于语音信号处理的通用工具.它定义了一个非常有用的时间和频率分布类, 其指定了任意信号随时间和频率变化的复数幅度. 实际上,计算短时傅里叶变换的过程是把一个较长的时间信号分成相同长度的更短的段, 在每个更短的段上计算傅里叶变换, 即傅里叶频谱.
短时傅里叶变换通常的数学定义如下:
其中,
DTFT (Decrete Time Fourier Transform) 为离散时间傅里叶变换.  其数学公式, 如下所示:
其中,  x(n) 为在采样数 n 处的信号幅度. ω~ 的定义如下:
实现时, 短时傅里叶变换被计算为一系列加窗数据帧的快速傅里叶变换 (Fast Fourier Transform, FFT)，其中窗口随时间 “滑动” (slide) 或“跳跃” (hop) 。
Python 实现
在程序中, frame_size 为将信号分为较短的帧的大小, 在语音处理中, 通常帧大小在 20ms 到 40ms 之间. 这里设置为 25ms, 即 frame_size = 0.025;
frame_stride 为相邻帧的滑动尺寸或跳跃尺寸, 通常帧的滑动尺寸在 10ms 到 20ms 之间, 这里设置为 10ms, 即 frame_stride = 0.01. 此时, 相邻帧的交叠大小为 15ms;
窗函数采用汉明窗函数 (Hamming Function) ;
在每一帧, 进行 512 点快速傅里叶变换, 即 NFFT = 512. 具体程序如下:
#-*- coding: utf8 -*-
importnumpy as npdef calc_stft(signal, sample_rate=16000, frame_size=0.025, frame_stride=0.01, winfunc=np.hamming, NFFT=512):#Calculate the number of frames from the signal
frame_length = frame_size *sample_rate
frame_step= frame_stride *sample_rate
signal_length=len(signal)
frame_length=int(round(frame_length))
frame_step=int(round(frame_step))
num_frames= 1 + int(np.ceil(float(np.abs(signal_length - frame_length)) /frame_step))#zero padding
pad_signal_length = num_frames * frame_step +frame_length
z= np.zeros((pad_signal_length -signal_length))#Pad signal to make sure that all frames have equal number of samples
#without truncating any samples from the original signal
pad_signal =np.append(signal, z)#Slice the signal into frames from indices
indices = np.tile(np.arange(0, frame_length), (num_frames, 1)) +\
np.tile(np.arange(0, num_frames* frame_step, frame_step), (frame_length, 1)).T
frames= pad_signal[indices.astype(np.int32, copy=False)]#Get windowed frames
frames *=winfunc(frame_length)#Compute the one-dimensional n-point discrete Fourier Transform(DFT) of
#a real-valued array by means of an efficient algorithm called Fast Fourier Transform (FFT)
mag_frames =np.absolute(np.fft.rfft(frames, NFFT))#Compute power spectrum
pow_frames = (1.0 / NFFT) * ((mag_frames) ** 2)returnpow_framesif __name__ == '__main__':importscipy.io.wavfileimportmatplotlib.pyplot as plt#Read wav file
#"OSR_us_000_0010_8k.wav" is downloaded from http://www.voiptroubleshooter.com/open_speech/american.html
sample_rate, signal = scipy.io.wavfile.read("OSR_us_000_0010_8k.wav")#Get speech data in the first 2 seconds
signal = signal[0:int(2. *sample_rate)]#Calculate the short time fourier transform
pow_spec =calc_stft(signal, sample_rate)
plt.imshow(pow_spec)
plt.tight_layout()
plt.show()
参考资料
1. DISCRETE TIME FOURIER TRANSFORM (DTFT). https://www.dsprelated.com/freebooks/mdft/Discrete_Time_Fourier_Transform.html
2. THE SHORT-TIME FOURIER TRANSFORM. https://www.dsprelated.com/freebooks/sasp/Short_Time_Fourier_Transform.html
3. Short-time Fourier transform. https://en.wikipedia.org/wiki/Short-time_Fourier_transform
4. Speech Processing for Machine Learning: Filter banks, Mel-Frequency Cepstral Coefficients (MFCCs) and What's In-Between. https://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html

起源

之前制作钢琴频谱识别时候，发现傅里叶变换的神奇之处，激起了内心的好奇。前几天在bilibili又发现了用傅里叶变换绘图的视频。决定自己做一个，顺便仔细学习一下傅里叶变换。

过程

从对傅里叶变换一窍不通，到反复学习推导过程，以及自己边做边发现问题，加上自己的思考和参考视频中的讲解，终于对傅里叶变换有了一个比较深刻的理解。

原理

对于一个波，对其进行傅里叶变换之后可以得到若干角速度、初相角、大小不同的复数向量。

复数向量可以以箭头形式在复平面内表示。将这些箭头首尾相连，其运动一端的轨迹便可以在复平面内形成二维图形。

而此二维图形与输入的波之间有一定联系。

依据此关系，最终可以用一定的波进行傅里叶变换之后，使变换结果在复平面内绘出对应于原波的图形。

成果

演示动图：









 

程序：

ps：

        1.需要注意的地方看程序注释。

        2.输入图片需要是单连通的，具体定义可以百度。对于一些多连通的图形，用ps或画图 稍作修改即可成为单连通图形（比如演示中的“A” 和音乐符号）。

        3.图形中尽量不要出现尖锐的角，由于寻找边界的算法很简陋（自己拍脑袋随便想了一下搞出来的），可能会寻到到尖锐处就无法继续。

 

function fft_draw()
    %picture=imread('C:\Users\zxy\Pictures\Matlab\doo.png');
    picture=imread('timg.png');
    %  1 ,  务必使用  白   底  黑  图  案  的  图片
    %  2，因为只能绘制外轮廓，所以 务必使用  单  连  通  图形
    %  3，图形务必居中（ 大致居中即可 ）。
    N_point=800;%最终绘制时的采样点数，越大越精细，一般400足够
    N_arrows=1;%显示的向量数；N_arrows=1或超过最大值时，自动取最大值
    show=1;% show=1时显示寻找边界的动画，为其他值时不显示
    cycle=2;% cycle=-1时无限循环显示动画，为n循环显示n次
    
    
    
    picture=RuiHua(picture);%锐化
    data_margin=ShiBie(picture, show);%识别边界
    Show(data_margin,N_point, N_arrows,cycle);

    function result=RuiHua(picture)%锐化，保证边界拾取的可靠性
        [width, high,~]=size(picture);
        for i=1:1:width
            for j=1:1:high
                if picture(i,j,1)<125
                    picture(i,j,:)=0;
                else
                    picture(i,j,:)=255;
                end
            end
        end
        result=picture;
    end

    function data=ShiBie(picture, show)
        %识别边界，得到边界点在图像中的坐标的数组（Xi,Yi）,其中i=1，2，3...n
        [width, high, ~]=size(picture);%获取图片宽度、高度
        high_find=round(high/2);%从图片1/2高度开始寻找边界点
        width_find=0;
        %边界识别的起始点
        for i=2:1:width-3
            if abs( picture(high_find,i,1)- picture(high_find,i+1,1))>125 && width_find==0
                %如果红色R（red）数组中，左侧点的亮度值 与 右侧点亮度值 的差 大于 125，则认为此处是边界。
                %并记录此处坐标（high_find，width_find）.记录的是边界上-内侧-的点的坐标
                width_find=i+1;
            end
        end
        margin=[high_find ,width_find];%边界点坐标数组，第一个点的坐标
        picture_show=imshow(picture);
        pause(2);
        %开始识别边界
        while picture(high_find,width_find ,3)-picture(high_find,width_find ,2)<100
            %在经锐化的黑白图中，每个点RGB值只可能是(255,255,255)或(0,0,0).
            %--如果出现例如(0,0,255)的情况，说明这个点已经被标记过。
            %--又因为图形区域是单连通区域。再次遇到已经被标记过的点，说明所有边界已经识别完成
            picture(high_find, width_find , 3)=255;
            %将B（RGB的R）数组值改为（255），RG保持不变（0，0），作为已经被记录的边界点的标志。
            if show==1
                picture(high_find-1:high_find+1,width_find-1:width_find+1 , 1)=255;
                %将识别点上下左右1的范围内的点的 R值全置为255，显示边界更明显
                set(picture_show,'Cdata',picture);
                pause(0.001);
            end
            for i=0:0.7854:5.5  % 0.7875 ~ pi/4
                x_find_1=width_find+round(cos(i ) );
                y_find_1=high_find+round( sin(i ) );
                x_find_2=width_find+round(cos(i + 0.7854) );
                y_fing_2=high_find+round( sin(i + 0.7854 ) );
                %围绕上一次的识别点，从右侧开始，逆时针方向依次对比相邻两个点的 B值
                %如果 B 值出现大幅下降，说明此处是边界，且较小值所在点是图形内侧。
                if picture(y_find_1, x_find_1, 2) -picture(y_fing_2, x_find_2, 2) >125
                    high_find=y_fing_2;
                    width_find=x_find_2;
                    margin=[margin ;high_find,width_find];%记录新的边界点的坐标
                end
            end
        end
        data=margin;
    end

    function Show(data_margin, N_point, N_arrows,cycle)
        [N,~]=size(data_margin);%边界点数
        data_point=zeros(N_point,2);
        
        arrows=plot([0 0],[0 0],'marker','.','color','b');hold on;%箭头
        orbit=plot([0 0],[0 0],'r');%轨迹
        
        xdata_X=zeros(1,2*N);
        ydata_X=xdata_X;
        xdata_Y=xdata_X;
        ydata_Y=xdata_X;
        
        data_margin=data_margin/100;
        data_margin(: ,1)=data_margin(: ,1) - data_margin( 1,1);%平移图形，使起/止点坐标为（0，0）
        data_margin(: ,2)=data_margin(: ,2) - data_margin( 1,2);
        data_X=fft(data_margin(:,1));%对 x 傅里叶变换
        data_Y=fft(data_margin(:,2));%对 y 傅里叶变换
        [n, ~]=size(data_X);
        ct=zeros(1,n);%各阶向量的初始角速度
        
        
        limt=max(abs(data_X)*0.7);%控制绘出图形的显示大小
        xlim([-limt limt]);
        ylim([-limt limt]*0.8);
        
        n=round(n/2);%阶数，傅里叶变换得到的是对称数组，取一半
        if N_arrows==1 || N_arrows>n
            N_arrows=n;
        end
        corner=zeros(n,2);
        
        while cycle==-1 || cycle>0
            if cycle>0
                cycle=cycle-1;
            end
            
            for n_point=1:1:N_point
                for i=1:1:n
                    ct(i)=ct(i)+(i-1)/N_point;%角速度
                    %理论上ct(i)= 0,pi, 2pi, 3pi ... npi
                    %在此比例为0：1：2：3：...n即可，具体大小只影响动画速度

                    if n_point==1 %--------------------------- x(t) ----------------------
                        if abs( data_X( i ) )>0
                            if imag(data_X( i) )>0
                                corner(i,1 )=acos( real(data_X(i))/abs(data_X(i))  );
                                %用傅里叶变换得到的复平面向量计算初始转角
                            else
                                corner( i ,1 )=-acos( real(data_X(i))/abs(data_X(i))  );
                            end
                        else
                            corner(i,1 )=0;
                        end

                        if abs( data_Y( i ) )>0 %   ---------------------- y(t) ----------------------
                            if imag(data_Y( i) )>0%初始转角
                                corner(i ,2 )=acos( real(data_Y(i))/abs(data_Y(i))  );
                            else
                                corner(i ,2 )=-acos( real(data_Y(i))/abs(data_Y(i))  );
                            end
                        else
                            corner(i ,2 )=0;
                        end
                    end


                    if i>1%  ------------- x(t) ---------------
                        %对于X而言，只需要x方向的运动，因此需要附加一个反向运动的向量 将y方向的运动抵消
                        xdata_X( i*2-1 )= 0.1*abs(data_X(i)) * cos( 2*pi*ct(i) + corner( i ,1 ) )+ xdata_X( i*2-2 );
                        xdata_X(  i*2 )=0.1*abs(data_X(i)) * cos( -2*pi*ct(i) - corner( i ,1 ) )+ xdata_X( i*2-1);
                        ydata_X( i*2-1 )=0.1*abs(data_X(i)) * sin( 2*pi*ct(i) +corner( i ,1 ) ) + ydata_X( i*2 -2 );
                        ydata_X( i*2 )=0.1*abs(data_X(i)) * sin( -2*pi*ct(i) -corner( i ,1 ) ) + ydata_X( i*2-1);
                    else
                        xdata_X(i)=0.2*abs(data_Y(i)) * cos( corner( i ,1 ) );
                        ydata_X(i)=0;
                    end

                    if i>1% ------------y(t) ----------------
                        %对于Y而言，只需要y方向的运动，因此需要附加一个反向运动的向量 将x方向的运动抵消
                        xdata_Y( i*2-1 )=-0.1*abs(data_Y(i)) * sin( 2*pi*ct(i) + corner( i ,2 ) )+ xdata_Y( i*2-2 );
                        xdata_Y(  i*2 )=-0.1*abs(data_Y(i)) * sin( -2*pi*ct(i) - corner( i ,2 ) )+ xdata_Y( i*2-1);
                        ydata_Y( i*2-1 )=0.1*abs(data_Y(i)) * cos( 2*pi*ct(i) +corner( i ,2 ) ) + ydata_Y( i*2-2 );
                        ydata_Y( i*2 )= 0.1*abs(data_Y(i)) * cos( -2*pi*ct(i) -corner( i ,2 ) ) + ydata_Y( i*2-1);
                    else
                        xdata_Y(i)=0;%
                        ydata_Y(i)=0.2*abs(data_Y(i)) * cos( corner( i ,2 ) );
                    end

                end
                x_draw=xdata_X(1:N_arrows) + xdata_Y(1:N_arrows) ;
                y_draw=ydata_X(1:N_arrows)+ydata_Y(1:N_arrows);
                data_point(n_point,:)=[y_draw(N_arrows) -x_draw(N_arrows)];
                set(arrows,'xdata',y_draw(2:N_arrows)  ,'ydata', -x_draw(2:N_arrows));
                set(orbit,'xdata', data_point(1:n_point,1) ,'ydata', data_point(1:n_point,2) );
                pause(0.009);
                drawnow;
            end
        end
    end
end