R - 訊號處理

# 訊息

## English

``````> install.packages('signal')
--- Please select a CRAN mirror for use in this session ---

Content type 'application/zip' length 255433 bytes (249 Kb)

package ‘signal’ successfully unpacked and MD5 sums checked

> bf <- butter(5, 0.2)

> library('signal')

Attaching package: ‘signal’

The following objects are masked from ‘package:stats’:

filter, poly

> bf <- butter(5, 0.2)
> freqz(bf\$b, bf\$a)
> freqz(bf)
> ch <- cheby2(5, 20, 0.2)
> freqz(ch, Fs = 100)
> zplane(ch)
> t <- seq(0, 1, by = 0.01)
> x <- sin(2*pi*t*2.3) + 0.25*rnorm(length(t))
> z <- filter(ch, x)
> plot(t, x, type = "l")
> lines(t, z, col = "red")
> grpdelay(ch, Fs = 100)
- Group delay (gd) calculated at 512 points.
- Frequencies (w) given in *Hz*.
gd          w
3.162407 0.00000000
3.162772 0.09765625
3.163869 0.19531250
3.165699 0.29296875
.......  .......

In grpdelay.default(filt\$b, filt\$a, ...) :
grpdelay: setting group delay to 0 at singularity
>```
```

``````## The R implementation of these routines can be called "matlab-style",
bf <- butter(5, 0.2)
freqz(bf\$b, bf\$a)
## or "R-style" as:
freqz(bf)
## make a Chebyshev type II filter:
ch <- cheby2(5, 20, 0.2)
freqz(ch, Fs = 100) # frequency plot for a sample rate = 100 Hz
zplane(ch) # look at the poles and zeros
## apply the filter to a signal
t <- seq(0, 1, by = 0.01) # 1 second sample, Fs = 100 Hz
x <- sin(2*pi*t*2.3) + 0.25*rnorm(length(t)) # 2.3 Hz sinusoid+noise
z <- filter(ch, x) # apply filter
plot(t, x, type = "l")
lines(t, z, col = "red")
# look at the group delay as a function of frequency
grpdelay(ch, Fs = 100)```
```
``````> library(signal)

Attaching package: ‘signal’

The following objects are masked from ‘package:stats’:

filter, poly

> data(package="signal")
> ls("package:signal")
[1] "an"            "Arma"          "as.Arma"       "as.Zpg"
[5] "bartlett"      "bilinear"      "blackman"      "boxcar"
[9] "butter"        "buttord"       "cheb1ord"      "chebwin"
[13] "cheby1"        "cheby2"        "chirp"         "conv"
[17] "decimate"      "ellip"         "ellipord"      "fftfilt"
[21] "FftFilter"     "filter"        "FilterOfOrder" "filtfilt"
[25] "fir1"          "fir2"          "flattopwin"    "freqs"
[29] "freqs_plot"    "freqz"         "freqz_plot"    "gausswin"
[33] "grpdelay"      "hamming"       "hanning"       "ifft"
[37] "impz"          "interp"        "interp1"       "kaiser"
[41] "kaiserord"     "levinson"      "Ma"            "medfilt1"
[45] "MedianFilter"  "pchip"         "poly"          "polyval"
[49] "remez"         "resample"      "roots"         "sftrans"
[53] "sgolay"        "sgolayfilt"    "specgram"      "spencer"
[57] "spencerFilter" "triang"        "unwrap"        "Zpg"
[61] "zplane"
> lsf.str("package:signal")
an : function (degrees)
Arma : function (b, a)
as.Arma : function (x, ...)
as.Zpg : function (x, ...)
bartlett : function (n)
bilinear : function (Sz, ...)
blackman : function (n)
boxcar : function (n)
butter : function (n, ...)
buttord : function (Wp, Ws, Rp, Rs)
cheb1ord : function (Wp, Ws, Rp, Rs)
chebwin : function (n, at)
cheby1 : function (n, ...)
cheby2 : function (n, ...)
chirp : function (t, f0 = 0, t1 = 1, f1 = 100, form = c("linear", "quadratic",
"logarithmic"), phase = 0)
conv : function (x, y)
decimate : function (x, q, n = if (ftype == "iir") 8 else 30, ftype = "iir")
ellip : function (n, ...)
ellipord : function (Wp, Ws, Rp, Rs)
fftfilt : function (b, x, n = NULL)
FftFilter : function (b, n)
filter : function (filt, ...)
FilterOfOrder : function (n, Wc, type, ...)
filtfilt : function (filt, ...)
fir1 : function (n, w, type = c("low", "high", "stop", "pass", "DC-0", "DC-1"),
window = hamming(n + 1), scale = TRUE)
fir2 : function (n, f, m, grid_n = 512, ramp_n = grid_n/20, window = hamming(n +
1))
flattopwin : function (n, sym = c("symmetric", "periodic"))
freqs : function (filt, ...)
freqs_plot : function (w, ...)
freqz : function (filt, ...)
freqz_plot : function (w, ...)
gausswin : function (n, w = 2.5)
grpdelay : function (filt, ...)
hamming : function (n)
hanning : function (n)
ifft : function (x)
impz : function (filt, ...)
interp : function (x, q, n = 4, Wc = 0.5)
interp1 : function (x, y, xi, method = c("linear", "nearest", "pchip", "cubic",
"spline"), extrap = NA, ...)
kaiser : function (n, beta)
kaiserord : function (f, m, dev, Fs = 2)
levinson : function (x, p = NULL)
Ma : function (b)
medfilt1 : function (x, n = 3, ...)
MedianFilter : function (n = 3)
pchip : function (x, y, xi = NULL)
poly : function (x)
polyval : function (coef, z)
remez : function (n, f, a, w = rep(1, length(f)/2), ftype = c("bandpass", "differentiator",
"hilbert"), density = 16)
resample : function (x, p, q = 1, d = 5)
roots : function (x, method = c("polyroot", "eigen"))
sftrans : function (Sz, ...)
sgolay : function (p, n, m = 0, ts = 1)
sgolayfilt : function (x, p = 3, n = p + 3 - p%%2, m = 0, ts = 1)
specgram : function (x, n = min(256, length(x)), Fs = 2, window = hanning(n), overlap = ceiling(length(window)/2))
spencer : function (x)
spencerFilter : function ()
triang : function (n)
unwrap : function (a, tol = pi, dim = 1)
Zpg : function (zero, pole, gain)
zplane : function (filt, ...)
> B1 <- c(0.2066,0.4131,0.2066)# Coefficients of numerator polynomial
> A1 <- c(1,-0.3695,0.1958)       # Coefficients of denominator polynomial
> H1z <- freqz(B1,A1,100)# Compute the transfer function
> # Filter H2(z)
> B2 <- c(0.894,-1.789,0.894)
> A2 <- c(1,-1.788,0.799)
> H2z <- freqz(B2,A2,100)
>
> # Filter H3(z)
> B3 <- c(0.42,0,-0.42)
> A3 <- c(1,-0.443,0.159)
> H3z <- freqz(B3,A3,100)
>
> # Filter h4(z)
> B4 <- c(0.5972,0.4425,0.5972)
> A4 <- c(1,0.4425,0.1584)
> H4z <- freqz(B4,A4,100)
> # Convenience function to draw multiple plots
> hPlot <- function(H){
+       text <- deparse(substitute(H))  # get the name of the filter for the title
+       c <- substr(text,4,4)
+       plot(H\$f,abs(H\$h),
+    col="red",
+                    ylim=c(0,1),
+            xlab="Normalized Frequency",
+            ylab="Magnitude",
+            main=paste("Filter H",c,"(z)",sep="")
+ )
+ }
>
> par(mfrow=c(2,2))
> plotList <- list(H1z,H2z,H3z,H4z)
> lapply(plotList,hPlot)
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

> #
> q()
> plot(H4z)# Look at the default plot
>```
```

page revision: 4, last edited: 25 Oct 2013 11:49