信号处理入门之旅

内容速览

readmatrix 函数从文件中导入数值数组

可以使用 readtimetable 函数将数据导入时间表.使用 "SampleRate" 选项计算每个采样的时间。

        mytbl = readtimetable("fn.txt","SampleRate",Fs) (Fs采样率是时间步长的倒数)

        "TimeStep" 时间步长  ,用 milliseconds、seconds 和 minutes 等函数向数值数组添加单位。 

可以使用 timerange 函数从时间表中提取一个范围：

lim = timerange(t1,t2)       "close"

获取某个时间范围内的所有数据，请使用冒号运算符：

tbl = tbl(lim,:)

查看信号的频率成分，可以使用 pspectrum 函数计算功率谱。

pspectrum(mysig)

为了放大有意义的地震活动，您可以在计算功率谱时设置频率范围。

pspectrum(mysig,"FrequencyLimits",[a b])

可以使用 pspectrum 函数获得功率谱和对应的频率。

[pwr,freq] = pspectrum(tbl)

如果 tbl 包含多个信号，则 tbl 是矩阵，其中每列对应一个信号的频谱。

开始创建分块布局之前，请定义布局中绘图的数量：

tiledlayout(2,3)  此代码创建一个 2 行 3 列的布局。

请使用 nexttile 函数填充该分块布局。

tiledlayout(2,3)
nexttile
plot(x,y)
nexttile
plot(x2,y2)

可以使用 synchronize 函数合并这些时间表。

t = synchronize(t1,t2)

地震信号均为均匀采样信号，因此在本练习中，您只需提供信号和重采样因子 p 和 q。 

y = resample(x,p,q)

p 和 q 是整数重采样因子。输出信号 y 具有 x 的 p/q 倍采样。

例如，

y = resample(x,2,3)

以原始采样率的 2/3 倍速率对 x进行重采样。

normalize 函数可用于对时间表进行归一化。

tbl = normalize(tbl)

默认情况下，normalize 函数使用 Z 分数方法。

要计算信号之间的互相关性，请使用 xcorr 函数。

[c,lags] = xcorr(x,y)

如果两个信号之间存在足够的相关性，可使用 finddelay 函数估计它们之间的相对延迟。

d = finddelay(x,y)

输入信号 x 和 y 应为数组，而不是时间表。输出 d 是采样个数。 

可以使用圆点表示法设置时间表属性：

tbl.Properties.PropName

将 harp 表的 StartTime 属性设置为 harpDelay。

harp.Properties.StartTime = harpDelay

通过使用 stackedplot 函数可以很方便地在一个时间表中显示多个变量。

stackedplot(tbl)

尝试在信号分析器中查看 quakesROI 信号：

signalAnalyzer(quakesROI)

semilogx 函数将创建一个绘图，其中 x 轴具有对数刻度。

semilogx(f,p)

可以直接使用 10*log10(p)，也可以使用 db 函数进行相同的计算。

db(p,"power")

semilogx(f,db(p,"power"))

一次只能创建一个信号的频谱图，但 quakes 表包含三个信号。您可以使用点索引分别输入信号和时间向量。

pspectrum(data.Signal,data.Time,"spectrogram")

可以使用 "MinThreshold" 选项去除低功率频率。

pspectrum(data,"spectrogram", ...
    "MinThreshold",n) 

可以用 cwt 函数创建尺度图。

cwt(sig,fs)

第一个输入是信号，第二个输入是采样率。

还可以使用 cwt 函数设置频率范围。

cwt(sig,fs, ...
    "FrequencyLimits",[a b])

cwt 函数没有 MinThreshold 选项，但您可以通过设置颜色图范围来实现相同的效果。

caxis([a b])

在给定频率上添加一条垂直线，有助于选择用于滤波的通带频率。

xline(pass)

可以使用 lowpass 函数进行低通滤波。

lowpass(tbl,pass)

要对时间表中的一个变量进行滤波，您可以使用冒号运算符 (:) 来获取给定变量名称的所有时间戳。

tbl(:,"VarName")

在 0.1 Hz 处对 quakes 中的 WANC 变量进行低通滤波。

lowpass(quakes(:,"WANC"),0.1)

要减少滤波信号中的高频成分，可以增大陡度。 

sig = lowpass(tbl,pass,"Steepness",s)

可以使用带通滤波器只保留该范围内的频率。

bandpass(tbl,[f1 f2])

在绘制低通和带通滤波信号之前，请将它们合并到一个时间表中。

tbl = timetable(t,sig1,sig2,...
    'VariableNames',["A" "B"])

plot后将 x 范围设置为从 2900 秒到 2950 秒

xlim([seconds(2900) seconds(2950)])

在基于频谱图进行任何计算之前，您首先需要获得数组（而不仅仅是图像）。

[p,f,t] = pspectrum(sig,time,"spectrogram")

三个输出是每段信号的频谱估计值、频率和时刻。

要找出哪些时间戳包含大量高频成分，您可以对每个时间戳的所有频率的功率求和。sum

psum = sum(p)

plot(t,psum)

为了突显峰值，您可以计算功率。

db(p,"power")

您可以使用“求局部极值”的交互式任务来自动识别信号中的局部最小值或最大值。