关于希尔伯特变换的 c语言实现

合集下载

emd 希尔伯特黄变换程序

emd 希尔伯特黄变换程序

(一)简单的EMD程序function imf = emd(x)% Empiricial Mode Decomposition (Hilbert-Huang Transform)% imf = emd(x)% Func : findpeaksx = transpose(x(:));%转置imf = [];while ~ismonotonic(x) %当x不是单调函数,分解终止条件x1 = x;sd = Inf;%均值%直到x1满足IMF条件,得c1while (sd > 0.1) | ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时,分量终止条件s1 = getspline(x1);%上包络线s2 = -getspline(-x1);%下包络线x2 = x1-(s1+s2)/2;%此处的x2为文章中的hsd = sum((x1-x2).^2)/sum(x1.^2);x1 = x2;endimf{end+1} = x1;x = x-x1;endimf{end+1} = x;% FUNCTIONSfunction u = ismonotonic(x)%u=0表示x不是单调函数,u=1表示x为单调的u1 = length(findpeaks(x))*length(findpeaks(-x));if u1 > 0, u = 0;else, u = 1; endfunction u = isimf(x)%u=0表示x不是固有模式函数,u=1表示x是固有模式函数N = length(x);u1 = sum(x(1:N-1).*x(2:N) < 0);u2 = length(findpeaks(x))+length(findpeaks(-x));if abs(u1-u2) > 1, u = 0;else, u = 1; endfunction s = getspline(x)%三次样条函数拟合成元数据包络线N = length(x);p = findpeaks(x);s = spline([0 p N+1],[0 x(p) 0],1:N);-------------------------------------------------------------------------------- function n = findpeaks(x)% Find peaks.找到极值% n = findpeaks(x)n = find(diff(diff(x) > 0) < 0);u = find(x(n+1) > x(n));n(u) = n(u)+1;------------------------------------------------------------------------------------------ ---------------------------------------------------------------------------------------- function plot_hht00(x,Ts)% 双边带调幅信号的EMD分解% Plot the HHT.% plot_hht(x,Ts)%% :: Syntax% The array x is the input signal and Ts is the sampling period.% Example on use: [x,Fs] = wavread('Hum.wav');% plot_hht(x(1:6000),1/Fs);% Func : emd% Get HHT.clear all;close all;Ts=0.0005;t=0:Ts:1; % 采样率2000HZ% 调幅信号x=sin(2*pi*t).*sin(40*pi*t);s1 = getspline(x);%上包络线s2 = -getspline(-x);%上包络线x1 = (s1+s2)/2;%此处的x2为文章中的hfigure;plot(t,x);xlabel('Time'), ylabel('Amplitude');title('双边带调幅信号');hold on;plot(t,s1,'-r');plot(t,s2,'-r');plot(t,x1,'g');imf = emd(x);for k = 1:length(imf)b(k) = sum(imf{k}.*imf{k});th = angle(hilbert(imf{k}));d{k} = diff(th)/Ts/(2*pi);end[u,v] = sort(-b);b = 1-b/max(b);% Set time-frequency plots.N = length(x);c = linspace(0,(N-2)*Ts,N-1);%figure;for k = v(1:2)plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3); hold on;set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 50]);%设置x、y轴句柄xlabel('Time'), ylabel('Frequency');title('原信号时频图');end% Set IMF plots.M = length(imf);N = length(x);c = linspace(0,(N-1)*Ts,N);for k1 = 0:4:M-1figurefor k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2});set(gca,'FontSize',8,'XLim',[0 c(end)]);title('EMD分解结果');endxlabel('Time');end(二)function imf = emd(x)% Empiricial Mode Decomposition (Hilbert-Huang Transform)% imf = emd(x)% Func : findpeaksx = transpose(x(:));%转置imf = [];while ~ismonotonic(x) %当x不是单调函数,分解终止条件x1 = x;sd = Inf;%均值%直到x1满足IMF条件,得c1while (sd > 0.1) | ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时,分量终止条件s1 = getspline(x1);%上包络线s2 = -getspline(-x1);%下包络线x2 = x1-(s1+s2)/2;%此处的x2为文章中的hsd = sum((x1-x2).^2)/sum(x1.^2);x1 = x2;endimf{end+1} = x1;x = x-x1;endimf{end+1} = x;% FUNCTIONSfunction u = ismonotonic(x)%u=0表示x不是单调函数,u=1表示x为单调的u1 = length(findpeaks(x))*length(findpeaks(-x));if u1 > 0, u = 0;else, u = 1; endfunction u = isimf(x)%u=0表示x不是固有模式函数,u=1表示x是固有模式函数N = length(x);u1 = sum(x(1:N-1).*x(2:N) < 0);u2 = length(findpeaks(x))+length(findpeaks(-x));if abs(u1-u2) > 1, u = 0;else, u = 1; endfunction s = getspline(x)%三次样条函数拟合成元数据包络线N = length(x);p = findpeaks(x);s = spline([0 p N+1],[0 x(p) 0],1:N);--------------------------------------------------------------------------------------------------------------------------------------------------------------- function n = findpeaks(x)% Find peaks.找到极值% n = findpeaks(x)n = find(diff(diff(x) > 0) < 0);u = find(x(n+1) > x(n));n(u) = n(u)+1;------------------------------------------------------------------------------------------ ---------------------------------------------------------------------------------------- function plot_hht00(x,Ts)% 双边带调幅信号的EMD分解% Plot the HHT.% plot_hht(x,Ts)%% :: Syntax% The array x is the input signal and Ts is the sampling period.% Example on use: [x,Fs] = wavread('Hum.wav');% plot_hht(x(1:6000),1/Fs);% Func : emd% Get HHT.clear all;close all;Ts=0.0005;t=0:Ts:1; % 采样率2000HZ% 调幅信号x=sin(2*pi*t).*sin(40*pi*t);s1 = getspline(x);%上包络线s2 = -getspline(-x);%上包络线x1 = (s1+s2)/2;%此处的x2为文章中的hfigure;plot(t,x);xlabel('Time'), ylabel('Amplitude');title('双边带调幅信号');hold on;plot(t,s1,'-r');plot(t,s2,'-r');plot(t,x1,'g');imf = emd(x);for k = 1:length(imf)b(k) = sum(imf{k}.*imf{k});th = angle(hilbert(imf{k}));d{k} = diff(th)/Ts/(2*pi);end[u,v] = sort(-b);b = 1-b/max(b);% Set time-frequency plots.N = length(x);c = linspace(0,(N-2)*Ts,N-1);%figure;for k = v(1:2)plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3); hold on;set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 50]);%设置x、y轴句柄xlabel('Time'), ylabel('Frequency');title('原信号时频图');end% Set IMF plots.M = length(imf);N = length(x);c = linspace(0,(N-1)*Ts,N);for k1 = 0:4:M-1figurefor k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2}); set(gca,'FontSize',8,'XLim',[0 c(end)]);title('EMD分解结果');endxlabel('Time');endplot_hht00(x,Ts)只画出了时频图,没有三维图。

希尔伯特(Hilbert)变换

希尔伯特(Hilbert)变换

希尔伯特(Hilbert)变换希尔伯特(Hilbert)变换是一种信号处理中常用的数学工具之一,主要用于将实数信号转化为复数信号,并提取出复信号的包络和瞬时相位等信息。

本文将对希尔伯特变换的基本概念、性质以及在信号处理中的应用进行介绍。

一、基本概念1. 复信号的生成在信号处理中,我们往往需要将一个实数信号变为一个复数信号,这可以通过对信号进行“解析”的方式来实现。

具体地,我们将实数信号x(t)通过一个信号处理器H(t)(即称为系统传递函数)得到一个复数信号X(t),即:X(t) = H(t) * x(t)其中,符号“*”表示对那些对应时间点处的信号进行点乘,即乘上相应的复数模长e^(jw),其中w为角频率,j为单位复数。

2. 复信号的包络和瞬时相位由于复数信号包含实部和虚部两个分量,其中实部和虚部分别表示原信号的信号值和90度相位移的信息。

因此,我们可以通过分别从复数信号中提取出它的实部和虚部,来获得原始信号的包络和瞬时相位两个信息。

具体的,假设我们有一个复数信号X(t) = x(t) + j*y(t),其中x(t)为实部,y(t)为虚部,则:信号的包络:A(t) = sqrt(x^2(t) + y^2(t))其中,atan2(y(t), x(t))表示y(t)/x(t)的反正切,但与通常的反正切最大的区别在于,它不仅考虑了y(t)/x(t)的值,而且也考虑了x(t)的符号,从而在所有象限范围内都具有唯一性。

3. 希尔伯特变换希尔伯特变换是一种用于从实数信号中构造复数信号的技术。

具体地,假设我们有一个实数信号x(t),那么它的希尔伯特变换y(t)定义如下:y(t) = H[x(t)] = P.\ I.C.\ \lim_{\varepsilon \to 0} \frac{1}{\pi}\int_{-\infty}^\infty \frac{x(t')}{t-t'-j\varepsilon} dt'其中,P和I.C.分别表示柯西主值和积分常数项。

Hilbert变换C源码

Hilbert变换C源码

Hilbert变换C源码/******************************************/ /* main.c/******************************************/#include <stdio.h>#include <math.h>#define N 500#define sq(X) ((X)*(X))int main(){int i;double x;double delta[N];double kappa[N];double y;double xmin = -10.;double xmax = 10.;double cd = -1.;double w = 2.;double h = (xmax - xmin) / N;for (i = 0; i < N; i++){x = 2. * (xmin + i * h - cd) / w;if (x < -1.)delta[i1] = 0.;else if (x < 1.)delta[i1] = sqrt(1. - sq(x));elsedelta[i1] = 0.;}hilbert(n, delta, kappa);for (i = 0; i < N; i++){x = 2. * (xmin + i * h - cd) / w;if (x < -1.)y = x + sqrt(sq(x) - 1.);else if (x < 1.)y = x;elsey = x - sqrt(sq(x) - 1.);(void) printf("%.3f %.3f %.3f %.3f\n", xmin + i * h, delta[i], k ppa[i], y);}return 0;}/******************************************//* hilbert.c/******************************************/void hilbert(int, double[], double[]);/******************************************//* hilbert.c/******************************************/#include <stdio.h>#include "hilbert.h"#define PI 3.14159265void hilbert(int n, double delta[], double kappa[]) {double d1, d2, d3, d4;int i1, i2;for (i1 = 0; i1 < n; i1++){kappa[i1] = 0.;for (i2 = 1; i2 < n; i2++){d1 = (i1+i2<N)? delta[i1+i2]: 0.;d2 = (i1-i2>=0)? delta[i1-i2]: 0.;d3 = (i1+i2+1<N)? delta[i1+i2+1]: 0.;d4 = (i1-i2-1>=0)? delta[i1-i2-1]: 0.;kappa[i1] -= 0.5 * (d1-d2) / i2 + 0.5 * (d3 - d4) / (i2 1);}kappa[i1] /= PI;}/***********************************************/* 说明/***********************************************The Hilbert transformThis package performs a numerical Hilbert transform. Download Compressed tar-fileorhilbert.c and hilbert.hHeader fileThe program must include the header file#include "hilbert.h"/***********************************************/* 说明kappa[i1] -= 0.5 * (d1-d2) / i2 + 0.5 * (d3 - d4) / (i2 1);}kappa[i1] /= PI;}}/***********************************************/* 说明/***********************************************The Hilbert transformThis package performs a numerical Hilbert transform. Download Compressed tar-filehilbert.c and hilbert.hHeader fileThe program must include the header file#include "hilbert.h"and then call hilbert() to perform the transform.UsageIf delta and kappa are arrays of n doubles, both arrays are al located by the maiprogram.The input data are in delta and the call hilbert(n, delta, kapp a) will return thHilbert transform of delta in kappa.n and delta are not modified.CompilationCompile the package using a C-compiler.If you use a makefile and this makefile looks like thisa.out: a.ob.oc.occ a.o b.o c.o...a.o: a.c d.hcc -c a.c...and you use the package in a.c, change the makefile to look like thisa.out: a.ob.oc.o hilbert.occ a.o b.o c.o hilbert.o...a.o: a.c d.h hilbert.hcc -c a.c...hilbert.o: hilbert.c hilbert.hcc -c hilbert.cDemomain.c uses the package to perform a Hilbert transform of a semi-ellipsis.Output from the program is a table.First column: xSecond column: The semi-ellipsisThird column: The Hilbert transform calculated by hilbert().Fourth column: The analytical Hilbert transform of the semi -ellipsis.。

希尔伯特黄变换获取时频谱, python

希尔伯特黄变换获取时频谱, python

希尔伯特黄变换(Hilbert-Huang Transform, HHT)是一种非线性、非平稳信号分析方法,能够有效地获取信号的时频谱信息。

在信号处理和振动分析领域,HHT被广泛应用于信号的时间-频率特征提取、故障诊断、模式识别等方面。

而Python作为一种功能强大的编程语言,为HHT的实现提供了便利条件。

下面将介绍希尔伯特黄变换的基本原理及其在Python中的实现。

1. 希尔伯特变换希尔伯特变换是对信号进行解析的一种数学方法,其核心是通过与原始信号相关的虚部信号来构建解析信号。

希尔伯特变换可以将实部信号与虚部信号相互转换,从而实现对信号的时域和频域分析。

希尔伯特变换的数学表示如下:\[H(x(t)) = P \left( \frac{1}{\pi t} \right) \ V \int_{-\infty}^{\infty} \frac{x(\tau)}{t-\tau} d\tau \]其中,\(x(t)\)为原始信号,\[H(x(t))\]为对应的希尔伯特变换,\(P\)表示柯西主值,\(V\)表示广义积分。

在时频分析中,希尔伯特变换可以用于提取信号的振幅和相位信息,从而实现时域和频域特征的全面分析。

2. 黄变换黄变换是由我国科学家黄次寅于1998年提出的一种基于希尔伯特变换的信号分析方法。

与传统的傅立叶变换和小波变换相比,黄变换更适用于非线性和非平稳信号的分析。

黄变换包括两个核心步骤:经验模态分解(EMD)和希尔伯特谱分析。

EMD是将复杂信号分解成若干个本征模态函数(EMD),而希尔伯特谱分析是在每个本征模态函数上进行希尔伯特变换,从而获取每个本征模态函数的时频特征。

3. 希尔伯特黄变换希尔伯特黄变换是将希尔伯特变换与黄变换相结合的一种信号分析方法。

希尔伯特黄变换主要包括以下步骤:1) 对原始信号进行EMD分解,得到若干个本征模态函数;2) 对每个本征模态函数进行希尔伯特变换,得到每个本征模态函数的时频谱信息;3) 将每个本征模态函数的时频谱信息相加,得到原始信号的时频谱分布。

c语言 希尔伯特曲线

c语言 希尔伯特曲线

c语言希尔伯特曲线
希尔伯特曲线是一种由德国数学家大卫·希尔伯特在1891年定
义的分形曲线,也被称为希尔伯特空间填充曲线。

它具有许多有趣的性质和应用,包括在图像压缩、信号处理和计算机图形学中的应用。

在C语言中,我们可以通过递归算法来实现希尔伯特曲线的绘制。

具体的实现过程如下:
1. 定义一个绘制函数,该函数接受四个参数:起始点的坐标、
线段的长度、绘制方向和递归的深度。

2. 在绘制函数中,首先根据当前的方向和线段长度计算出终止
点的坐标。

3. 如果递归深度为0,表示已经到达最底层,直接绘制一条线
段连接起始点和终止点。

4. 如果递归深度不为0,那么需要先递归绘制下一层的曲线。

具体的方法是将线段长度除以2,然后将当前方向逆时针旋转90度
分别得到三个新的方向,接着分别以起始点和终止点为起点,递归调用绘制函数。

5. 最后,在所有递归调用完成后,再绘制一条线段连接起始点
和终止点。

通过以上的实现方法,我们就可以在C语言中实现希尔伯特曲线的绘制了。

这不仅可以帮助我们深入理解分形图形的生成原理,还可以展示C语言的递归和图形绘制能力。

- 1 -。

证明希尔伯特变换是正交的

证明希尔伯特变换是正交的

证明希尔伯特变换是正交的希尔伯特变换的正交性在信号处理中有着非常重要的应用。

通过对其进行研究和证明,我们可以更好地了解希尔伯特变换的性质以及在信号处理中的实际应用。

下面,我们将为大家介绍证明希尔伯特变换是正交的详细过程。

首先,我们需要明确什么是希尔伯特变换。

希尔伯特变换是一种线性运算,用于将一个信号转换为其复解析信号,即在傅里叶变换的基础上,将实部和虚部的信号分别反转并取相反数,用于分析信号在时域和频域中的行为。

那么,如何证明希尔伯特变换是正交的呢?步骤一:定义内积首先,我们需要定义内积。

内积是在函数空间中定义的一种运算,用于度量两个函数相似程度的大小。

对于两个实函数f(x)和g(x),其内积定义如下:( f , g )= ∫f(x)g(x)dx对于两个复函数,其内积定义如下:(f , g)= ∫f*(x)g(x)dx其中,f*(x)是f(x)的共轭复数。

步骤二:证明正交性根据内积的定义,我们可以证明希尔伯特变换是正交的。

首先,我们需要证明希尔伯特变换是有限范围内的,即其信号在无穷远处趋于0。

根据奇偶性,我们可以证明实部和虚部对应的傅里叶变换值分别关于频率轴对称,因此,它们的线性组合即希尔伯特变换的傅里叶变换值也是关于频率轴对称的,这就保证了其在无穷远处趋于0。

接下来,我们需要证明希尔伯特变换的内积等于零。

我们可以将希尔伯特变换的实部和虚部表示为:h(x)=f(x)cos(wx)-g(x)sin(wx)h*(x)=f(x)cos(wx)+g(x)sin(wx)其中,f(x)和g(x)是两个实函数,w是一个常数。

我们可以将希尔伯特变换的内积表示为:(h , h*)= ∫[f(x)cos(wx)-g(x)sin(wx)][f(x)cos(wx)+g(x)sin(wx)]dx将其展开,得到:(h , h*)= ∫[f(x)^2+g(x)^2]cos^2(wx)dx +∫[f(x)^2+g(x)^2]sin^2(wx)dx根据三角函数的性质,cos^2(wx)+sin^2(wx)=1,因此,上式可以简化为:(h , h*)= ∫[f(x)^2+g(x)^2]dx由于f(x)和g(x)都是实函数,因此其和的平方和是非负的,即(h , h*)≥ 0。

1希尔伯特变换的基本原理

1希尔伯特变换的基本原理

1希尔伯特变换的基本原理希尔伯特变换(Hilbert transform)是一种非常重要的信号处理技术,它在时间域和频率域之间建立了一种特殊的变换关系,可以通过提取信号的相位信息来分析信号的时频特性。

本文将详细介绍希尔伯特变换的基本原理。

一、定义与表达式希尔伯特变换首先由德国数学家大卫·希尔伯特(David Hilbert)提出,他建立了一个衍生(Analytic)函数的概念。

对于一个实值信号函数x(t),它的希尔伯特变换H{x(t)}可以表示为:H{x(t)} = \frac{1}{\pi} \int_{-\infty}^{\infty}\frac{x(\tau)}{t-\tau} d\tau其中,H{x(t)}是实值信号的希尔伯特变换,x(t)是原始信号,t是时间变量。

希尔伯特变换可以通过对信号的频谱进行处理实现,首先对原始信号进行傅里叶变换得到频谱X(f),然后将频谱进行处理后再进行逆傅里叶变换得到希尔伯特变换。

具体来说,对于一个实值信号x(t),它的傅里叶变换为X(f),那么它的希尔伯特变换H{x(t)}可以表示为:H{x(t)} = IFT \{ -j \cdot sign(f) \cdot X(f) \}其中,IFT 表示逆傅里叶变换,sign(f)是频率变量的符号函数。

二、频谱分析希尔伯特变换的一个重要应用是信号的频谱分析,通过分析信号的相位信息来了解信号的时频特性。

希尔伯特变换可以提取信号的边带频率信息,从而反映信号的局部属性。

对于一个实值信号x(t),它的频谱X(f)可以分解为实部和虚部:X(f) = X_r(f) + j \cdot X_i(f)其中,X_r(f)和X_i(f)分别是实部和虚部的频谱函数。

希尔伯特变换可以通过将频谱的虚部部分置零来获得信号的解析信号。

解析信号是一种由实信号和其希尔伯特变换构成的复信号表示,它具有可分辨信号的相位信息的特点。

三、希尔伯特变换的性质希尔伯特变换具有许多重要的性质,其中最重要的性质是希尔伯特变换的平移性质和相位信息的提取。

希尔伯特变换 时域做法

希尔伯特变换 时域做法

希尔伯特变换在时域的实现方法可以通过以下步骤进行:
1. 确定输入信号f(t)。

2. 构造一个冲击响应为1/πt的单位冲激响应函数g(t)。

3. 将输入信号f(t)与单位冲激响应函数g(t)进行卷积运算,得到输出信号h(t)。

这个过程在时域中相当于将输
入信号通过一个滤波器,该滤波器的冲击响应为g(t)。

4. 输出信号h(t)即为输入信号f(t)的希尔伯特变换结果。

需要注意的是,上述方法中的单位冲激响应函数g(t)是连续时间希尔伯特变换的情况。

对于离散时间希尔伯特变换,单位冲激响应函数会有所不同,但基本的实现思路是相似的。

另外,在实际应用中,由于计算机只能处理离散信号,因此需要对连续信号进行采样和离散化处理。

在离散化后,可以利用离散时间傅里叶变换(DTFT)或离散余弦变换(DCT)等方法来计算希尔伯特变换。

总之,希尔伯特变换在时域的实现方法主要是通过卷积运算将输入信号与一个特定的单位冲激响应函数进行滤波处理,从而得到输出信号。

这个过程可以看作是输入信号通过一个特定的滤波器或系统。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

关于希尔伯特变换的c语言实现
最近毕设需要用到希尔伯特变换的知识,今天做完之后决定还是记录下思路。

当然是数字信号的希尔伯特变换
上面是连续信号的希尔伯特变换,离散的应该也能根据上面写(没现成的图片,懒得编辑公式了)。

这里打算采用使用卷积的方法来计算。

由于希尔伯特变换的传输函数的傅里叶变换是H(w)= -j w&gt;=0
j w&lt;0
所以我们可以先求原始信号的离散傅里叶变换E(K),然后按照下面的公式就可以求出希尔伯特变换之后的信号的离散傅里叶变换Eh(K)。

然后对Eh(K)求反傅里叶变换就可得到我们需要的信号的希尔伯特变换信号。

下面贴代码思路
先建立一个复信号的结构体:
typedef struct {
Float64 r; /* 实部*/
Float64 i; /* 虚部*/
} CPX;接着是离散傅里叶变换的函数第一个参数dir代表正变换和反变换。

void DFT(int dir,int framelen,CPX *signal,CPX *dft_s)
{
int i,k;
double arg;
double cosarg,sinarg;for(i=0;i&lt;framelen;i++)
{
arg=-dir*2.0*3.141592654*(double)i/(double)framelen; for(k=0;k&lt;framelen;k++)
{
cosarg=cos(k*arg);
sinarg=sin(k*arg);
dft_s[i].r+=(signal[k].r*cosarg-signal[k].i*sinarg);
dft_s[i].i+=(signal[k].r*sinarg+signal[k].i*cosarg); }
}
/*返回数据*/
if(dir==-1)
{
for(i=0;i&lt;framelen;i++)
{
dft_s[i].r=dft_s[i].r/(double)framelen;
dft_s[i].i=dft_s[i].i/(double)framelen;
}}
}上主函数
int main(void)
{
CPX *signal;
CPX *dft_s;
CPX *hdft_s;
CPX *hsignal;signal=calloc(framelen,sizeof(CPX)); //
原始信号
dft_s=calloc(framelen,sizeof(CPX)); // 原始信号的傅里叶变换
hdft_s=calloc(framelen,sizeof(CPX)); // 希尔伯特变换的离散傅里叶变换
hsignal=calloc(framelen,sizeof(CPX)); // 希尔伯特变换后信号DFT(1,framelen,signal, dft_s); //求原始信号傅里叶变换for(i=0;i&lt;framelen;i++) //求出希尔伯特变换信号的傅里叶变换
{
if(i&lt;=framelen/2)
{
hdft_s[i].r=dft_s[i].i;
hdft_s[i].i=-dft_s[i].r;
}
else
{
hdft_s[i].r=-dft_s[i].i;
hdft_s[i].i=dft_s[i].r;
}
}
DFT(-1,framelen, hdft_s,hsignal); //利用反傅里叶变换求出希尔伯特变换信号}
最后这个hsignal就是我们要的希尔伯特变换信号了。

相关文档
最新文档