Cellular communication uplink I can see, from my apartment, a cell tower.

From what I have learned, in a time division duplex system the uplink and downlink communication are kept separate in time. Uplink from some time, then downlink, then uplink, and so on. I have also learned that during uplink all of the mobiles transmit simultaneously to the base station.

So this leads to my question:

in a TDD system, how do all of the mobiles keep track of time so that they are all in synch?

For example, from what I know, my mobile and my neighbors mobile both are transmitting to the base station during the uplink time slot, but how do they keep their timing in synch with each other?

BlackMath hit the nail on its proverbial head:

  Using control channels from the base station.

Typically, the whole circus works like this:

Before a user equipment can register at a base station, it has to look for that cell's "beacon". Whenever it sees that, it waits a time $T_\text{wait}$ after the beacon, at which new devices are allowed to say "hey, I'm here". Because these devices can't know of each other yet, that phase isn't really TDD per se (it's more like ALOHA, with the random element being delay do to distance); especially, because these devices are at different distances from the base station (which they can't know already), it's not time-slotted. For that reason, the "hello" phase is not as tightly timed as regular TDD slots. That means it "wastes" a lot of uplink time, and for that reason, the network operator might define these phases to happen relatively rarely.

Assuming the new user equipment successfully registered within that "hello" phase (i.e. didn't collide with another new device), the base station now knows that there's a new device, and this devices cumulative delay $T_\text{beacon} + T_\text{distance} + T_\text{wait} + T_\text{distance}$. Since to the base station, $T_\text{beacon}$ and $T_\text{wait}$ are known, it can easily deduct $T_\text{distance}$. It then uses the downlink to tell the new device that, "hey, you need to be $T_\text{timing advance}$ earlier than you are, to fit into our TDD scheme. Your timeslots to send data are 1, 3, 22 and 128. Welcome to my cell!".

The user equipment adjusts its internal TDD timekeeping using that timing advance. During operation, the cell tower continuously informs the user equipment whether it needs to do adjustments of that timing advance, which would, for example, be the case if the user moved away or towards the cell tower.


Bilinear transformation of continuous time state space system I'm trying to understand the derivation of the bilinear transform for a set of continuous time state-space matrices. I've found plenty of websites which list steps to perform the conversion (here 1 or here 2 or here 3) - but haven't found any derivation or insight into how they came about. When I attempt to do it myself, I keep getting terms in z which I can't seem to get rid of.

To make the problem clear, if I define the bilinear transform as:

$$ s=\frac{\alpha\left(z-1\right)}{z+1} $$

Where alpha is usually given as 2/T.

The transfer function of a continuous-time state-space system can be given as:

$$ H(s) = \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} + \mathbf{D} $$

How do I derive the expressions for the matrices of the discrete time version bilinear transformed model:

$$ H(z) = \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d} $$

When I make the variable change, I keep getting extra terms in z which I don't know what to do with. I'm kinda hoping this is just me failing to have a strong grasp of matrix operations. Any insight (including if I am going about this completely the wrong way) would be really appreciated!

--- Edit with answer.

Thanks to Klaz for posting an answer to this first. I managed to work through the problem today and thought I would post my own alternate working here.

Given the typical state-space equations:

$$ s\mathbf{Q}(s) = \mathbf{A}\mathbf{Q}(s) + \mathbf{B}\mathbf{X}(s) $$
$$ \mathbf{Y}(s) = \mathbf{C}\mathbf{Q}(s) + \mathbf{D}\mathbf{X}(s) $$

Focus on the state update equation first and make the bilinear substitution:

$$ \frac{\alpha\left(z-1\right)}{z+1} \mathbf{Q}(z) = \mathbf{A}\mathbf{Q}(z) + \mathbf{B}\mathbf{X}(z) $$

We manipulate this expression to make the LHS $z\mathbf{Q}(z)$:

$$ \alpha\left(z-1\right) \mathbf{Q}(z) = \mathbf{A}\left(z+1\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

$$ z\left(\alpha\mathbf{I}-A\right)\mathbf{Q}(z) = \left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

$$ z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

This instantly gives us $\mathbf{A_d}$ as:

$$ \mathbf{A_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right) $$

We know that the linear $(z+1)$ term of $\mathbf{X}(z)$ can be applied wherever $\mathbf{Q}(z)$ is used. We apply this in the output expression to give the new output discrete equation as:

$$ z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\mathbf{X}(z) $$
$$ \mathbf{Y}(z) = \mathbf{C}\left(z+1\right)\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$
$$ \mathbf{Y}(z) = \mathbf{C}z\mathbf{Q}(z) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

We can now pull out $\mathbf{B_d}$ as:

$$ \mathbf{B_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B} $$

Because we know an expression for $z\mathbf{Q}(z)$, we can formulate $\mathbf{Y}(z)$ as:

$$ \mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}\mathbf{Q}(z) + \mathbf{B_d} \mathbf{X}(z) \right) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

$$ \mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right)\mathbf{Q}(z) + \left(\mathbf{C}\mathbf{B_d} + \mathbf{D}\right)\mathbf{X}(z) $$

Giving:

$$ \mathbf{C_d} = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right) $$
$$ \mathbf{D_d} = \mathbf{C}\mathbf{B_d} + \mathbf{D} $$

I've had the same question last week, but I've managed to find how to derive it (getting rid of those $z$ terms is indeed tricky). I will give here detailed demonstration of how to arrive to the result given in 1 (with, in your notation, $\alpha = 2 \lambda$).

So we define our new discrete-time function transfer as

$$
\begin{array}{rcl}
H_d(z) &=& H(\frac{\alpha\left(z-1\right)}{z+1}) \\
&=& \mathbf{C} \left[ \frac{\alpha\left(z-1\right)}{z+1}\mathbf{I}-\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\alpha} \mathbf{C} \left[ \left(z-1\right)\mathbf{I}- \frac{z+1}{\alpha}\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \left(\mathbf{I} - \frac{1}{\alpha} \mathbf{A} \right) - \left(\mathbf{I} + \frac{1}{\alpha} \mathbf{A}\right) \right]^{-1}\mathbf{B} + \mathbf{D} \\
\end{array}
$$

For sake of notation, let $\mathbf{P} = \mathbf{I} - \frac{1}{\alpha} \mathbf{A}$ and $\mathbf{Q} = \mathbf{I} + \frac{1}{\alpha} \mathbf{A}$ (important remark useful later on: $\mathbf{P} + \mathbf{Q} = 2 \mathbf{I}$). Then

$$
\begin{array}{rcl}
H_d(z) &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{P} - \mathbf{Q} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \mathbf{P}^{-1} \mathbf{B} + \mathbf{D} \\
&=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \left( \sqrt{\frac{2}{\alpha}}  \mathbf{P}^{-1} \mathbf{B} \right) + \mathbf{D} \\
&=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D}
\end{array}
$$
where we defined
$ \mathbf{A_d} = \mathbf{P}^{-1} \mathbf{Q} $
and
$ \mathbf{B_d} = \sqrt{\frac{2}{\alpha}}  \mathbf{P}^{-1} \mathbf{B} $.
Furthermore, we have

$$
\begin{array}{rcl}
\frac{z+1}{\sqrt{2 \alpha}} \mathbf{C}
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z + 1 \right) \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z \mathbf{I} + \mathbf{I} + \mathbf{A_d} - \mathbf{A_d} \right) \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \left( \mathbf{I} + \mathbf{A_d} \right) \right] \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \mathbf{P}^{-1} \left( \mathbf{P} + \mathbf{Q} \right) \right] \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right]
\end{array}
$$

which gives us

$$
\begin{array}{rcl}
H_d(z)
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right] \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\
&=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} + \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\
&=& \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d}
\end{array}
$$
with
$ \mathbf{B_d} = \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} $
and
$ \mathbf{D_d} = D + \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} $.
Replacing $\mathbf{P}$ and $\mathbf{P}$ gives the results in 1.

PS: the same idea can be used to prove the discrete-time state-space representation found using  Generalized Bilinear Transform ($s \leftarrow \alpha \frac{z - 1}{\beta z + \left(1-\beta\right)}$) or First-order Holder methods.


Why do we always use square kernels for filters? In image processing, why do we always use odd square kernels(3x3, 5x5, 7x7, etc.) for filters?

Why can't we use kernels of size 3x1, 5x3, etc. that are rectangular kernels? Also, why do we not prefer even kernels (2x2, 4x4, 6x6, etc.)?

If the kernel is odd then you can center it on top of a pixel, which is nice. Lots of kernels are symmetric, so that's probably why they are often square. You can certainly use non-square or even kernels though.


Converting from a phase response to a group delay response and vice-versa I'm trying to figure out how to convert any phase response to the corresponding group delay response.

Is there a way to do that, and vice-versa if possible?

The group delay of a filter is defined as minus the change in the phase response with respect to frequency.  If the phase response of a filter is $\Phi(\omega)$, the corresponding group delay $\tau_g$ is given by:

\begin{equation}
\tau_g = -\frac{d\Phi(\omega)}{d\omega}
\end{equation}

In Matlab code, the group delay of a 4th order Butterworth filter can be calculated like so:

[b a] = butter( 4, 0.25 ); % design the filter
[Hz fVec] = freqz( b, a ); % compute the frequency response and angular frequencies
phi  = unwrap(angle(Hz)); % extract phase response
taug = -diff(unwrap(angle(Hz)))/((fVec(2))); % compute group delay
figure; plot( fVec(2:end), taug );  % plot the result

The phase response must be unwrapped to avoid $\pi$ to $-\pi$ discontinuities.  Because the difference between adjacent frequency bins is the same, we need only scale the -diff(phi) by 1/fVec(2).  Finally, because the diff command produces a vector one less than its input, the group delay values are plotted against fVec(2:end) and not fVec.

Alternatively, you could use the matlab command grpdelay, though you're likely to learn more from the previous!

figure; grpdelay(b,a);


Remove noise from linear signal This signal represents the distance in an experimental device. As the speed is constant, the signal should be linear (untill it drops). However, for some reason the device showed interferences as shown in the figure below:

The question is how can I remove this high amplitude interferences, making the signal as linear as possible?

I've already tried the butterworth filter and read something about the Chebyshev, but I don't think it applies.

Any suggestions?

There are a few simple things that you could try:

sliding-window median filtering with an appropriate window size
linear-regression in a sliding window
robust linear regression to the total signal, if you know the signal duration etc well enough

EDIT: Further idea after a comment from OP on another answer

Since you want the slope at 160000 be the same as before (i.e. no flat section in your signal), you could use a robust linear regression technique such as RANSAC with an appropriate number of sample points to analyze/filter your data.


Help me with BM3D algorithm's aggregation scheme BM3D  is considered state of the art in Image Denoising. I am trying to implement the algorithm in Matlab. Right now, I am done with the 3D wavelet transform in the first step, and to be precise, the equation 3 in the aforementioned paper. What I want to do now is to carry out the equation 5 on my de-noised patches in 3D domain. Can someone please explain the equation 5, the aggregation scheme, a bit more lightly?

BM3D intuition

The intuition behind BM3D is the same as other non-local algorithms: an image is made by the juxtaposition of many, self-repetitive small image elements called patches. Thus, when looking at small enough pieces, it is possible to find enough similar patches that can be combined and denoised together.

The different non-local algorithms are different in the way used to combine these similar patches.
BM3D relies on sparsity. This principle has become widely accepted now, and basically states that for many useful signals (including images) there exists at least one frame where this signal can be decomposed into a small number of meaningful components.

For images, DCT and wavelet basis are obvious candidates for the sparse frame: JPEG / JPEG2k compression rely on the transformation of the input images into 8-by-8 blocks followed by retaining the strongest coefficients, which is indeed a sparsifying transform.

When noise is added to a signal, its frequencies get mixed or added to the frequencies of the noise-free signal. For example, white noise adds non-zero values to the whole spectrum. Thus, finding a sparse approximate of a noisy signal (in a good basis of course) yields usually a good denoising performance.

BM3D mechanics

The 3D transform of BM3D is actually a 2D+1 transform.

Algorithmically, this is what happens (I'll describe the mathematical intuition afterwards):

make a stack of similar patches. All the patches of the stack are supposed to be small variants of the same ideal (noise free) patch. The variant term is important, you'll see why later;
take a 1D transform (DCT / wavelet) along one of the horizontal axes;
take a 1D transform (DCT / wavelet) along the other. What you have now is a stack of transformed patches;
now, take the transform of choice (DCT / wavelet / Wiener filtering) along the vertical axis. What you are doing here is taking a transform (DCT / wavelet) or directly computing a denoised value (Wiener filtering) for each point in the vertical column;
(optional) if you used DCT / wavelet transform, apply a sparsifying transform, usually hard thresholding.
Then, compute the denied patches by inverting the different transforms (in reverse order);
now, back project the result to the original patch locations. Patches can (and will) overlap when you denies the whole image. This raises an important practical question: how do you fuse the results on overlapping areas? Which denoised patch is the most reliable? BM3D uses Eq. 5 to fuse the results by a weighted averaging: the weight of a back-projected patch is inversely proportional to the variance in the patch stack. When the variance of the stack is high, it means that the patches did not really come from the same ideal one, and thus the result is marked as less reliable.

Mathematically, the intuition is that most patches should be sparse in some good domain (typically, DCT domain for $8 \times 8$ patches).
BM3D also supposes that patches in a stack are variants of a same patch, and thus the variations between these transformed patches and the ideal patch (in the transformed domain) should also be sparse. This assumption is enforced by the 3rd transform and the hard thresholding step.


SIFT - Taylor Expansion I'm trying to implement SIFT algorithm.
After I found the min/max points in DOG images I need to find the real min/max subpixels by taylor exapnsion:

$$
\begin{equation}
D(\mathbf{x}) = D + \frac{\partial D^T}{\partial \mathbf{x}}\mathbf{x} + \frac{1}{2}\mathbf{x}^T\frac{\partial^2D}{\partial \mathbf{x}^2}\mathbf{x}
\end{equation}
$$

http://en.wikipedia.org/wiki/Scale-invariant_feature_transform

I'm not sure how I do it.
I understand the idea behind using taylor expansion.

Consider one-dimensional function $f(x)$.

The first order taylor expansion is $f(x_0+h) \approx f(x_0) + f'(x_0)h$

The second order taylor exapnsion is $f(x_0+h) \approx f(x_0) + f'(x_0)h + \frac 1 2 f''(x_0)h^2$

Now we expand three-dimensional function.
$$ D(\mathbf{x_0}+\mathbf h) \approx D(\mathbf{X_0}) + \bigg(\frac{ \partial D}{\partial \mathbf x}\bigg)^T\bigg|_{\mathbf x=\mathbf x_0}\mathbf h + \frac 1 2 \mathbf h^TH(\mathbf x)\mathbf h $$
, where
$$
\frac {\partial D}{\partial \mathbf x} = \begin{bmatrix} \frac {\partial D}{\partial x} \\ \frac {\partial D}{\partial y} \\ \frac {\partial D}{\partial \sigma}\end{bmatrix} =
\begin{bmatrix} \frac {D(x+1,y,\sigma) - D(x-1,y,\sigma)}{2} \\ \frac {D(x,y+1,\sigma)-D(x,y-1,\sigma}{2} \\ \frac {D(x,y,\sigma+1) - D(x,y,\sigma-1)}{2} \end{bmatrix}
$$
$$ H(\mathbf x) = \begin{bmatrix} D_{xx} & D_{xy} & D_{x\sigma} \\
D_{yx} & D_{yy} & D_{y\sigma} \\
D_{\sigma x} & D_{\sigma y} & D_{\sigma\sigma} \end{bmatrix} $$

Additionally, $D_{xy}$
$$= \frac { \frac{D(x+1,y+1,s)-D(x-1,y+1,s)}{2} - \frac{D(x+1,y-1,s)-D(x-1,y-1,s)}{2} } 2 $$

Let's go back.
$$D' = \bigg ( \frac{ \partial D}{ \partial \mathbf x} \bigg ) ^T + H(\mathbf x) \mathbf h = 0 $$ $$\mathbf h = - H^{-1}(\mathbf x) \bigg ( \frac{ \partial D }{\partial \mathbf x} \bigg ) ^T $$

Assume that $\mathbf{x_0}$ is discrete feature point. Iteratively we can get $\mathbf h$.

If $\mathbf{x_0}$ is the closest point to $\mathbf h$, then we select $\mathbf{x_0}+\mathbf h$. (Imagine that 3-dimensional space)

My english ability is poor and there may be some errors.


Gradient echo MRI. Do the gradients have to be the same duration? I currently got hold of the book "Handbook of MRI Pulse Sequences" and am looking into computing the duration, gradient strength and ramp times for the phase encoding gradient, readout prehasing and the slice prephasing gradients.

Now, every GE pulse sequence diagram that I look at, shows these gradients to be of the same duration i.e. they start and end at the same time during the sequence. Now, going through some of the calculations in the book for a scanner with some maximum gradient strength, slew rate, receiver BW, field of view and the k-space matrix size, I can compute the optimal gradient strengths and the durations for these gradients.

However, the optimal gradients are not the same duration. Does this mean I have made a mistake in my calculations or do they overlap these gradients in the pulse sequence diagrams for convenience?

The question is what optimal is. It is correct that very often certain gradients do have the same duration in sequence diagrams - and actual imaging sequences (e.g. slice rewinder, read-out prewinder, and phase-encoding gradient). Normally, you would compute the longest duration for these gradients and recalculate the others to be of the same length by stretching them in time and reducing their amplitude. This may have benefits, e.g. the peripheral nerve stimulation may be reduced, and also the noise that the MRI creates. Also it reduces the stress for the gradient coils, i.e. heating and vibration. Certainly it is not always the best idea to stretch all gradients: When the fall below a certain amplitude, there may be inaccuracies due to system imperfections. Also the parasitic diffusion weighting of the gradients may be an issue in certain applications.

So the main reasons for certain design choices are:

reduction of patient discomfort
reduction of stress for the MR-system
reduction of image artifacts


methods of crest factor minimization i am studying the methodes CF minimization. and  have a doubt.

in each methods there are steps with DFT and IDFT or like this FFT and IFFT.
Everyone provides the fact, that they use DFT and IDFT for calculation muti-sin with new phase and then take DFT or FFT to present in time domain.

But what is the purpose of dft and idft in minimization of the CF?
Are their function  to represent the signal in F -or T- domain?

Typically the goal of crest factor minimization is to find the time signal with the small-ish crest factor for a given amplitude spectrum in the frequency domain. If your goal is different, please update your question.

The only "dial" you have is the phase spectrum. Hence the primary output of a CF minimization process is a phase spectrum. With the exception of a few special cases there are no good analytical solutions for this (at least as far as I know), so most methods are iterative or heuristic.

If you choose a random phase, you typically end up with a crest factor of about 4 to 5, which is simply a consequence of the central limit theorem. Your amplitude distribution will be roughly Gaussian, which in theory has infinite crest factor. However, in practice, amplitudes higher than 5 are extremely unlikely (around 1e-6) so they rarely occur or the error resulting by simply clipping them off of negligible

Iterative methods work the following way:

Initialize with a random phase
Go to time domain
Clip time domain signal slightly below the target crest factor
Go back to the frequency domain and repair the amplitude spectrum
Repeat until you reached your target crest factor or you time out

The idea here, is that every iteration step drives the phase spectrum closer to one which represents a lower crest factor. This works quite well with reasonably "broad" spectra and you may be able to get to about 1.2 or so.

Sparsely populated spectra are more tricky since there are less degrees of freedom.

An alternative method is to construct the time domain signal directly by concatenating chunks of sine waves where the dwell time at each frequency is proportional to the target energy at this frequency. Then repair any spectral errors in the frequency domain. This will create more "sweep like" signals and tends to works better for spectrally sparse signals. CFs will come about 1.6 or so. In the special case of a white signal this is equivalent to the "squared phase" approach


Determine the type of filter from the transfer function i got that system $x_n \to x_n-x_{n-1}$, so $h_n=[.....,0,1,-1,0,...]$, with $h_0=1$ and $h_1=-1$, so the transfer function given by:

$$\sum_{i=-\infty}^{\infty} h_ne^{-jwn} = h_0e^{-jw(0)}+h_1e^{-jw(1)}=1-e^{-jw}$$. How to knwo the type of filter, if it is HPF, LPF, etc.

I appreciate your help.

That is an FIR filter transfer function.

Plot the frequency response to determine if it is HP, BP, LP.

(Hint: compare the magnitude at w=0 and w=pi to get a quick idea)

This seems like a homework problem, so I don't think it is appropriate to supply a complete answer.


Higher order filters I know how higher order filter affects the system, theoretically. The slope of frequency Response becomes steeper as its order increases.  I had some speed signals from work. I wanted to investigate more about this, So i designed a $N^th$ order $PT1^n$ filter(multiple PT1 filter concatenated in parallel) and a $N^th$ order Butterworth filter. I just cannot understand the stuff I read theoretically applies here.

1st to 3rd order $PT1^n$ filter:

1st to 3rd order butterworth filter:

These are real signals. It seems to me that as the order is being increased, it's increasing the phase shift to increase the noise reduction. I want to figure out if increasing the order of the filter is worth the hassle, as the order increases the complexity increases and the needs more computational power too.

Main purpose of my filtering is to smoothen the signal, remove those steps. I want to know how muchbenifit will i get by increasing the order of the filter, mathematically and practically

Main purpose of my filtering is to smoothen the signal, remove those steps. I want to know how much benifit will i get by increasing the order of the filter, mathematically and practically

There are basically two parameters for filters like a Butterworth lowpass: Its order and the cutoff frequency. Both have influence on the response.

Setting the cutoff frequency lower will make the response more smooth. Setting it higher will make the response follow more closely the steps because more of the frequencies are retained and less are attenuated.

The filter order controls how selective the filter is / how "sharply" it separates the passpand (the frequencies that are preserved) from the stop band (the frequencies that are rejected. Inbetween the passband and the stopband you have the transition band which will get narrower the higher order you use. In the time domain higher order filters tend to make the response more wavy. The filter's impulse response oscillates more and its envelope decays more slowly. Choosing the filter order is kind of a trade-off between frequency selectivity and time selectivity. By time selectivity I mean that a tiny pop in your source signal will affect a larger window of time in the filtered signal.

Without knowing exactly what you are shooting for, it's hard to suggest anything. I hope this explanation of how the two tuning knobs affect the response will enable you to move into the right direction.

If you want to do noise reduction it's generally good to know what the PSD of your signal and the noise part looks like. In your case, the noise is highly correlated and depends too much on the signal. You basically have truncation without any prior dithering which would otherwise decouple the noise error from the signal. So, in your case it's hard to predict the PSD of the error. But given both spectra, Wiener filtering could get you closer to the solution. Then, a good filter's cutoff frequency would depend on the point where the signal PSD crosses the noise's PSD. The order would depend on the difference of the PSD slopes between signal and noise in a loglog domain. So, for example, given a brownish signal mixed with white noise, the difference of the PSD slopes in the loglog domain would be 2. Divide this by two (--> 1) and use that as a filter order. The most appropriate cutoff frequency setting would be to make the -6 dB point the frequency where the signal PSD curve intersects the noise PSD curve.

If you care about the phase response, you could apply the IIR filter bidirectionally. Once in the forward direction and then in the backward direction. This cancels the phase response and squares the magnitude response. To compensate for this you'd reduce the filter order by about a factor of two. This might be a good approach if you do this kind of filtering offline. The phase and group delay of higher order IIR filters is higher like you observed.


Fourier series in continous time domain I want to ask Question about the Fourier series in continuous time domain.
I am following signal and systems 2nd Edition by Alan Oppenheim.
I have confusion in understanding the statement that
Specifically, suppose that $x(t)$ is real and can be represent in the form 3.25. then since $x^*(t) = x(t)$, we obtain

$$x(t) = \sum^{+\infty}_{k\ =\ -\infty} a^*_k e^{-jk\omega_0t}$$

Then it means the equation 3.25 is for both Real and imaginary?
Equation 3.25

$$x(t) = \sum^{+\infty}_{k\ =\ -\infty} a_k e^{jk\omega_0t}$$

Under certain conditions, a $T$-periodic function can be represented by its Fourier series

$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{jk\omega_0t}\tag{1}$$

with $\omega_0=2\pi/T$. The function $x(t)$ can be complex-valued, i.e. have non-zero real and imaginary parts. Note that generally the Fourier coefficients $a_k$ are also complex-valued (even for real-valued $x(t)$).

Now if $x(t)$ is real-valued, i.e., $x(t)=x^*(t)$, you get

$$x(t)=x^*(t)=\sum_{k=-\infty}^{\infty}a^*_ke^{-jk\omega_0t}=
\sum_{k=-\infty}^{\infty}a^*_{-k}e^{jk\omega_0t}\tag{2}$$

Comparing $(2)$ with $(1)$ you see that for real-valued $x(t)$ you get the condition $a_k=a^*_{-k}$, i.e. the Fourier coefficients show Hermitian symmetry. This condition is not satisfied for general complex-valued functions $x(t)$. However, the representation $(1)$ is valid no matter if $x(t)$ is real-valued or complex-valued.


Analyzing audio by looking at its time-frequency points in real- imaginary plane I am performing Short-time Fourier transform on an audio signal (I dont know how to attach the audio file here). The sound was recorded in anechoic chamber and the source is supposed to come from left. Lets say the result of STFT has 'N' number of time frames and 'nfft' FFT coefficients. My goal is to look at the plurality of time-frequency points in a certain frequency band over a plurality of time frames, graphically represented in real and imaginary parts in the complex plane. This is the code:

[x, fs] = audioread('booth.wav');   % load a .wav file x1 = x(:,1);
                                    %left ear

                          % bandpass filter to work on good frequencies only:
fc=[300 500]; %because my signal's frequency is in this range
[b,a] =butter(4,fc/(fs/2),'bandpass'); % Butterworth filter of order 4

x1 = filter(b,a,x1);

                                      %perform STFT:
x=x1(:);  %make sure it's column
xlen=length(x); %length of signal
wlen=2^nextpow2(0.04*fs);  %window length
hop=(3/4)*wlen;       %hop size
numFrames=1 + fix((xlen - wlen)/hop); %number of frames
nfft=2^nextpow2(wlen);

if mod(nfft,2)==0        %even
  N=nfft/2+1;
elseif mod(nfft,2)==1    %odd
  N=(nfft+1)/2;
end

window=hamming(wlen,'periodic'); %hamming window
Sl = zeros(numFrames,N); %preallocate stft for left ear
indx=0;

for i=1:numFrames
  xw=x(indx+1:indx+wlen).*window;   %multiplying signal to window
  X=(fft(xw,nfft))./wlen;  %fourier transform
  Sl(i,:)=X(1:N);    %keep half of it
  indx=indx+hop;   %hop and go to next time frame
end

f=(0:N-1)*fs/nfft;  %frequency vector
                    %correcting DC and Nyquist
if rem(nfft, 2)
  Sl(2:end, :) = Sl(2:end, :).*2;
else
  Sl(2:end-1, :) = Sl(2:end-1, :).*2;
end

                                %Now plotting the points:
figure(1);
hold on;
grid on;
for j=150:447   %I just want to look at these frequency indices
  plot(real(Sl(:,j)),imag(Sl(:,j)),'Ob');
  title(['F=', num2str(f(j)), ' hz'])
  xlabel('real'); ylabel('imaginary');
  axis([-1 1 -1 1]);
  plot(-1:1,[0,0,0],'k', 'LineWidth', 1);
  plot([0,0,0],-1:1,'k','LineWidth',1)
  pause;
end

If I look at the points at a frequency that I know my fft should have value at it, the points all fall on zero in the complex plane (if not on zero,very close to it). Why?
I tried doing RMS normalization, but the points now are not in the [-1 1] range in complex plane and they are all over the place! I am not sure if I can do something like peak normalization and Rms normalization together?

This should be a comment but comments have limittaitns

I don't kno if this fixes your problem because I'm not sure what you are doing but using the handel file that comes with matlab

Modifying your code

figure(1);
ii=intersect(find(f>300 ),find(f<500 ));
for j=ii %I just want to look at these frequency indices
plot(real(Sl(:,j)),imag(Sl(:,j)),'Ob'); %grid on;

title(['F=', num2str(f(j)), ' hz']);

produced something.

Hope this helped


Invariance of curvature in Curvature Scale Space PDF is here

I have read this PDF. But i don't know what the Invariance means in page 3 and 5. Calculation/estimation of κ is NOT invariant and Invariant to viewpoint. What's the meaning of these sayings? When i read "Invariant to viewpoint", i think it means the k value is same when rotation and translation. But after i read "Calculation/estimation of κ is NOT invariant" ,i have no idea. Then i try the code what i got here,to caculate the k value.

Code is here

And i have found that, the k value is different when shape have rotated which have the same sigma. Then i changed the sigma, k value is also different while the shape is same. When i set the sigma as 0.6667, the curvature of straight line is 1.6 not 0, why? The curvature of straight line should be zero,right? I don't why the value is 1.6. Then I tried to rotate the straight line, the curvature changed!!new value is -6.6. The papge said to find zero crossing to get concave-convex. But here curvature of straight line is not 0 is 1.6? Can somebody help me. Thank you very much!

OK, so let's look at the curvature of some standard objects.

First, a circle:
$$
x(s) = \cos(s)\ ;\ \dot{x}(s) = -\sin(s)\ ;\ \ddot{x}(s) = -\cos(s)\\
y(s) = \sin(s)\ ;\ \dot{y}(s) = \cos(s)\ ;\ \ddot{y}(s) = -\sin(s)\\
$$
so the curvature is:
$$
\kappa(s) = \frac{\sin(s)\sin(s) + \cos(s)\cos(s)}{(\sin^2(s) + \cos^2(s))^{3/2}} = 1
$$

Rotating the circle gives:
$$
x'(s) = \cos(\theta) x(s) - \sin(\theta) y(s)\\
y'(s) = \sin(\theta) x(s) + \cos(\theta) y(s)\\
$$
so that
$$
\dot{x}'(s) = -\cos(\theta)\sin(s) -\sin(\theta)\cos(s) = -\sin(\theta + s)\\ \ddot{x}'(s) = -\cos(\theta)\cos(s) +\sin(\theta)\sin(s) = -\cos(\theta + s)\\
\dot{y}'(s) = \sin(\theta)\sin(s) +\cos(\theta)\cos(s) = \cos(\theta - s)\\
\ddot{y}'(s) = \sin(\theta)\cos(s) - \cos(\theta)\sin(s) = \sin(\theta - s)\\
$$
Even thought it looks like I got a sign wrong somewhere, this will also evaluate to 1.

So, at least for a circle, the formula is rotationally invariant.

The trouble comes when you try to implement this.

function [kappa,xdot,ydot,xdotdot,ydotdot] = curvature(x,y,delta_s)
    xdot = diff([x x([1 2])])/delta_s;
    xdotdot = diff([xdot])/delta_s;
    ydot = diff([y y([1 2])])/delta_s;
    ydotdot = diff([ydot])/delta_s;

    xdot = xdot(1:length(xdot)-1);
    ydot = ydot(1:length(ydot)-1);

    kappa = (xdot.*ydotdot - xdotdot.*ydot)./sqrt(xdot.^2 + ydot.^2 + 10^9).^3;

endfunction

The only way to get the "derivatives" of your parameterized curves will be to take differences.  For the circle, this isn't too bad, but for more complex shapes, rotations might make these values a little different --- especially if the shapes are generated from two different sets of measurements.

For example the plot below shows:

The original curvature (note problem at end of range)
The curvature of a circle rotated
The curvature of a noisy circle
The curvature of a rotated noisy circle

// 25885
delta_s = .01;
s = -%pi:delta_s:(%pi-delta_s);
N = length(s);

x = cos(s);
y = sin(s);

xn = cos(s) + 0.001*rand(1,N,'norm');
yn = sin(s) + 0.001*rand(1,N,'norm');

[k,xd,xdd,yd,ydd] = curvature(x,y,delta_s);
[kn,xdn,xddn,ydn,yddn] = curvature(xn,yn,delta_s);

theta = %pi/4;
rot = [cos(theta), -sin(theta); sin(theta) cos(theta)];

all = rot*[x;y];
xdash = all(1,:);
ydash = all(2,:);

alln = rot*[xn;yn];
xdashn = alln(1,:);
ydashn = alln(2,:);

[k2,xd2,xdd2,yd2,ydd2] = curvature(xdash,ydash,delta_s);
[k2n,xd2n,xdd2n,yd2n,ydd2n] = curvature(xdashn,ydashn,delta_s);

clf
subplot(221)
plot(s,k)
title('Curvature of a noiseless circle')
subplot(222)
plot(s,k2)
title('Curvature of a ROTATED noiseless circle')
subplot(223)
plot(s,kn)
title('Curvature of a NOISY circle')
subplot(224)
plot(s,k2n)
title('Curvature of a ROTATED NOISY circle')


Compensating for group delay in FIR with symmetrical coefficients Suppose I have a sinusoidal signal of 50Hz sampled at 6.4kHz and I wished to filter it at:

Band pass: 160Hz, Bands stop: 800Hz, Apass:0.1db, Astop: 106db

For a FIR filter order of 40 I would get 41 number of taps. The group delay would thus be (41-1)/2 = 20.

Am i wrong in thinking that if I extended (/zero padded at the end of) the data block size to be filtered by the estimated group delay above, I would be compensating for the group delay?

For example if I had 128 samples (corresponds to a full sine cycle at 50Hz-Fs=6.4Khz) to be filtered and added the group delay 128+20=148. Now, if I were to filter these 148 samples with the above stated filter, would I not get a signal with no group delay?

Thanks
Alex

No.

There is no way you can compensate for a  broad-band group delay without violating causality. If you add zeros at the end of your signal, the group delay stays constant and if you add zeros in front of it, the group delay goes up.

At some frequencies it may look that there is no group delay (or phase shift), but that's just an illusion: the phase of the filter at this frequency just happens to be an integer multiple of two pi. If you look at any other frequency, it won't be true any more.


Can an EEG be re-sampled to fix a poor voltage resolution? I’ve made quite a bad mistake in ca. 14 EEG recordings – I recorded at 10uV resolution @ 5000 Hz instead of 0.1 uV @ 500 Hz. I’m conducting an ERP experiment, and the signals of interest are on the order of ~5 uV. I’m wondering if there is any way to up-sample the voltage given I have way more time series data points than I need..? Some sort of interpolation?

I’ve seen a number of posts on up-sampling from say 500 hz to 1000 hz, but not sure if the principal is the same?

Is there a way to trade off the sample rate against the poor resolution? Maybe by means of interpolation, or an averaging method to apply directly to the raw data**

This represents about 42 hours of recording time and I’m anxious to know if I can recover any usable data from these recordings, or if I have to try to get participants back in (there is a treatment enrollment deadline which means I can’t simply acquire more data).

The bit depth basically indicates how accurately your analog-to-digital converter records and reproduces the signal (Fig. 1). A higher bit depth means that more subtle fluctuations in the waveform are more faithfully reproduced in the digitized signal (source: Presonus)

As the answer from Bryan alludes to, data cannot be generated de novo, unfortunately. Data can be processed, yes, but data acquired at a certain resolution cannot be enhanced. Of course, as Strongbad indicates, you can generate an indefinite number of steps by upsampling, that's not an issue, but that will simply leave you with an estimate of possible intervening values; the real values are simply not there and are lost. Look at Fig. 1. If we would upsample the data, the points between successive samples would be interpolated, seemingly increasing your voltage steps, but the resulting trace will not be any better than it was before (only bigger :-).

Note that the resolution in the time-domain (sample rate) will not help you in any way. The amplitude (y) step size in Fig. 1 will not change, even if you has used an infinitely high sample rate.

Fig. 1. Illustration of bit depth. source: Presonus


Relationship between Bandwidth and Bit Rate In one of my classes I have learned that if we look at transmitting a binary signal with ASK, we can get the bandwidth (i.e. the width of the main lobe in the frequency domain) with BW = 2 * Fb, where Fb is the bit-rate (i.e. 1/T, where T = 'length' of one bit).

I am a little confused as to how this BW <-> Fb relation comes to be. Intuitively I understand that if we assume a carrier frequency of 2GHz and transmit a signal with a bitrate of 5 bits/s, we simply modulate the amplitude of our 2GHz carrier signal, without any frequency changes. Why, then, do we have a lobe in the frequency domain with a width that is greater than 1? As I understand it, there is only one frequency, the 2GHz signal?.. I know my understanding is wrong, but intuitively I do not understand how we have a variance in frequency, independently of what our bitrate is. How is modulating the amplitude of a 2GHz signal 5 times a second different from doing it 20 times a second when viewed in the frequency domain?

Assume you have a pulse of width $T$ which is your bit.

$$w(t)=1\quad for\quad-\frac T2\le t\le\frac T2\quad and\; 0\;elsewhere$$
 Your bit in frequency domain is:

$$W(f)=Sinc(T.f)$$

Look up Sinc function and you can see $Sinc(T.f)$ main lobe is between $\pm \frac 1T$. As your bit is widened in time, it occupies less bandwidth in frequency domain($T$ increases, $\frac 1T$ decreases) and vice versa. This is a very important result to understand and remember for the rest of your life. It's the reasoning behind relationship of frequency $BW$ and pulse duration $T$. Carrier frequency has nothing to do with it, it just shifts your entire signal spectrum up in frequency to make it suitable for transmission in the physical medium(channel) of your choice.

You should review Fourier transform and its properties.  Then look up Convolution theorem, which tells if you multiply signals in time, they get convolved in frequency domain. Then learn about Sinc function and its relationship to Rect function.


Find power of a sum of sinusoids We got this question to solve:

Calculate the power of the signal:
$$s(t) = 8\cos\left(20\pi t-\frac \pi4\right) + 4\sin(15\pi t)$$

Now, I thought of two approaches :

Use Parseval theorem, so first find Fourier transform, then calculate the required integral, but then we get a impulse function in it and I don't know how to integrate square of $\delta$.

Using the normal formula to calculate the power of a signal, and again in that, we get 3 terms, $\sin^2$ which is easy, then $\cos^2$ again simple, and then $\sin\cdot \cos$, which I know how to solve, but will be lengthy and I think there should be some another proper way to solve this question.

Any guidance would be appreciated.

HINT:
For real-valued periodic signal $x$ with period $T_0$ Parseval gives you
$$
P=\frac 1T_0\int_{-T_0/2}^{T_0/2} x^2(t)\mathrm{d}t = \sum_{n=-\infty}^{\infty}\lvert c_n\rvert^2
$$
Where $c_n$'s are Fourier coefficients. You can then use the properties in Equations $(1)$ and $(2)$.
$$
\mathcal F \big\{A\cos\left(2\pi f_0 t + \phi\right)\big\} = \frac A2\big[e^{j\phi}\delta\left(f - f_0\right) + e^{-j\phi}\delta\left(f + f_0\right)\big]\tag{1}
$$
$$
\mathcal F \big\{A\sin\left(2\pi f_0 t + \phi\right)\big\} = j\frac A2\big[e^{-j\phi}\delta\left(f + f_0\right)-e^{j\phi}\delta\left(f - f_0\right)\big]\tag{2}
$$
EDIT:
For Equation $(1)$ for instance, applying Parseval above you have the following:
\begin{align}
P&= \sum_{n=-\infty}^{\infty}\lvert c_n\rvert^2\\
&= \lvert c_{-1}\rvert^2 +  \lvert c_{1}\rvert^2\\
& = \frac{A^2}{4}\cdot 1  + \frac{A^2}{4}\cdot 1\tag{$\lvert e^{\pm j\phi} \rvert = 1$} \\
\end{align}


Does FFT magnitude of a sound drop to 1/10 of its value if sound level is decreased by 20dB? I understand decibel = 20 * log (magnitude/reference value). To truly appreciate the meaning of this, if I have a single frequency tone that corresponds to a certain FFT magnitude shown on a PC that's constantly recording. If I lower the tone volume by 20 dB (say, by monitoring with a SPL meter), then I should see the FFT magnitude that's 1/10 of its previous value right?

I'm only trying to understand this from a theoretical perspective, ignoring complicating real life factors.

Thanks

-20dB corresponds to $1/10$ of the amplitude or $1/100$ of the energy. That's the same thing, since energy (all else being equal) is proportional to the square of the amplitude.

It also corresponds very roughly to about $1/4$ perceived loudness but that's an entirely different can of worms


decorrelating detector In this paper decorrelating detector we have the decorrelating detector. With $K$ users, the correlation matrix of the signature waveform is $R=S^TS$ with $S=[s_1,...,s_K]$. Moreover $size(R)=K \times K$

This is my question: In the article (p $679$), we have: $d_1=\sum_{i=1}^K [R^{-1}]_{1,i}s_i$. But, how can we make this product, the elements do not have the same dimensions? In fact, $size(R^{-1})=K \times K$ and $size(s_i)=N \times 1$ with $N \neq K$

It all works out because $[\mathbf{R}^{-1}]_{1, i}$ is just a number. It is the number in row $1$ column $i$ of the matrix $\mathbf{R}^{-1}$. So the product $[\mathbf{R}^{-1}]_{1, i}\mathbf{s}_i$ returns a vector of length $N \times 1$. See the line after Equation 10.


How to calculate transfer rate and needed bandwidth I have two questions:

1- If I want to transfer a photo via bluetooth, how should I find how long it takes to transfer? I mean the formula between the image size, bluetooth bandwidth and time.

2- If I want to transfer video in real time, how should I calculate the bandwidth is needed to transfer the video via bluetooth/WiFi? What's the formula between video resolution(VGA,2K,4K,etc) and transferring bandwidth if we want to do it in real-time?

For example: If I have 1920*1080*30(fps) video, I think I should multiply it by 3 or 4(respect to colored image type) then I have about 180~240 million bytes to transfer per second. So it means I need about 200MBps bit-rate! And I only should use 802.11n protocol for that!

Did I calculate the things right?

EDIT: I found this calculator that I am not sure but I think it's the thing I was looking for (https://www.cctvcalculator.net/en/calculations/bandwidth-calculator/). what do you think about this? Why does it calculate local bandwidth bigger than network bandwidth? it's strange for me!

Your calculations are wrong. You have misunderstood the concept. The files may contain an amount of bits as you calculated but transferring that information is a bit different.

First of all, you send information as PACKETS. Here is an example packet type which is being used in Bluetooth.

The only place you send the actual information is called PAYLOAD. As you can see, payload is occupying a smaller area than you expected in a packet. The size of payload can smaller or larger depending which type of packet you are using, or which type of standards you are working with etc. And your bit rate is determined by which modulation type you use, how many bits your one symbol contains etc.

And not to forget, if your conditions are not perfect (which is the general situation) you WILL send some packets which is corrupt and you need to send them again. Which also reduces your actual bit-rate.

So there is different calculations for different combinations of the different techniques.

I hope the answer is clear.

As always,
Have a nice day.


Determining the pattern direction of a bitmap image How can I obtain average direction of the pattern shown in the figure below. It is the direction of the red arrow relative to the yellow (horizontal) line. Any ideas for an approach? I couldn't figure out a way to approach. This is a spatio-temporal image created from a video. Thank you.

Solution is using structure tensors

from skimage.feature import structure_tensor

Axx, Axy, Ayy = structure_tensor(gray_image, sigma=1)
tan2t = (2*Axy)/(Axx-Ayy)
t=np.arctan(tan2t)/2
data = t.ravel()

# find the mode of `data` array (can use histogram). it is the direction.

another question and answers from other community


What happens when we reconstruct a signal after its Fourier transform is done and some frequencies are altered? Since we loose time information when we take a Fourier transform, what happens if we alter a few frequencies from the transform and then reconstruct the signal, do we get the altered frequency at the exact place the original frequency was(in time). Would it be beneficial if we use Wavelet transform or FFT? If yes, how?

The spatial content is carried in the phase of the Fourier transform and not in magnitude. As Marcus pointed out, the difficulty with Fourier transform is that it doesn't allow us to immediately analyze the spatial information while still looking at the frequency. There are methods such as windowed FT which try to make that a bit better.

However, if you alter frequencies in the magnitude and reconstruct, you would get filtering effects. For instance, if we cut the low frequencies, than it would be a high-pass filter, and likewise, cutting high bands would cause a low-pass filter effect. This is analogous to applying, e.g. Mean/Gaussian filter to an image.

Regarding Wavelets, they try to address this time/frequency issue I mentioned above. Wavelets are limited in time and frequency. We slide the wavelets in time domain to construct the Wavelet domain. This provides resolution in the time domain.

On the application side, Wavelet denoising is shown to be superior to Fourier domain denoising. But, such superiority is application dependent, as for some applications one doesn't need that spatial analysis.

You could find some more information on the differentiation here.

Here is an example, based on this one, where we do a low-pass filter by keeping only the low frequency values (both positive and negative):

% make our noisy function
t = linspace(1,5,1024);
x = -(t-2).^2  + 2;
y = awgn(x,0.5);
F = fft(y,1024);

freqRange = 20; % range of frequencies we want to preserve

rectangle = zeros(size(F));
rectangle(1:freqRange+1) = 1;       % preserve low positive frequencies
y_half = ifft(F.*rectangle,1024);   % +ve low-pass filtered signal
rectangle(end-freqRange+1:end) = 1; % preserve low negative frequencies
y_rect = ifft(F.*rectangle,1024);   % full low-pass filtered signal

figure;
plot(t,y,'g--');
hold on, plot(t,x,'k','LineWidth',2);
hold on, plot(t,y_half,'b','LineWidth',2);
hold on, plot(t,y_rect,'r','LineWidth',2);
legend('noisy signal','true signal','positive low-pass','full low-pass','Location','southwest');

The resulting signals look like this:

You could see here that the altered frequencies effect the whole signal - wherever a high frequency existed gets smoothed. The spatial content is preserved in reconstructing the signal, because FT is a loss-less transformation, but we didn't get the chance to modify the exact spatial locations - the whole signal is effected by a change to the FT magnitude.


How to quantify area of interest in an image? How to find the percentage of a certain colour in a image accurately? For example  I am actually trying to find the percentage of barnacles in this image over the red part of the hull.
http://thumb1.shutterstock.com/display_pic_with_logo/1498544/150482510/stock-photo-barnacles-barnacles-perched-on-the-red-buoy-150482510.jpg

Given that the size of the barnacles is not expected to be considerably larger than the area of the hull being imaged at a certain shot, you can simply deduce the actual barnacle coverage area by counting how many pixels belong to the hull and then multiplying by a "pixel size" that you can derive by including an item of known length in the image.

The key assumption behind this (simple) technique is that the barnacles and the hull exist on the same plane (As if the barnacle layer was a "sticker" on the hull).

Obviously, this is not the case because the barnacles are three dimensional and "stick out" of the hull. For the same reason, the angle at which a certain part of the hull is imaged might also introduce an additional error factor.

You can control for all this by imaging relatively large areas of the hull (i.e. the camera is at a relatively long distance from the hull), taking care for the camera to be pointing vertically to as much of the hull as it is possible (this of course conflicts with the "take a picture at a long distance" and some experimentation will be required). Finally, you will also need to have good control of the lighting of the shot and for this reason, you might have to use a flash bulb (even in day-light) to flood the scene uniformly with light and avoid shadows.

Having taken the image, you can then convert its colour space from Red, Green, Blue (RGB) to Hue, Saturation, Brightness (HSB) and examine the Hue component's histogram for the uniform colour of the hull (here, Red and because of this, relatively easy to pick out). Again, this assumes that the hull is painted using just one colour.

Here is the process in Gimp:

The three small B&W images at the botoom of the screenshot are the HSB components respectively. The histogram at the top right is derived from the H component (the first B&W image at the bottom left-hand side of the screenshot).

As you can see, there is a clear peak towards the beginning of the H spectrum corresponding to the red colour (Please see the H band in this image to understand how the horizontal axis of the histogram now "corresponds" to hue value).

If there was an object of known length on the hull (an iron waterline marking relief for exaple or a pen or a coin) we could derive the pixel size and through this, the actual surface of the whole image. By then counting how many pixels "belong" to the peak we could deduce the actual surface occupied by the barnacles.

For example:

Image Dimensions: $450 \cdot 229 = 103050$ pixels

Pixels corresponding to red: $8494$

(Please note:You probably cannot see this because of the size of the image but once you have your histogram in GIMP (Colours->Info->Histogram) you can mark an area on the histogram and read how many pixels it describes. This is how this number was derived.).

Barnacle coverage area (%): $\frac{103050-8494}{103050} * 100 \approx 91.7\% $

All you need now is an object of known size to convert this to $cm^2$

In terms of automating this process, it can be done using a number of techniques and platforms (e.g. Python, Octave, etc) but for a more accurate answer reagrding those, a little bit more detail would be required regarding the application.

Hope this helps.


Why is C++ so popular for Image Processing? I have been a student of Image Processing for 3 years now and have been using C++ as the programming language to implement different Image processing operations. But I don't really why in particular C++ is so popular for Image processing rather than other languages like Java and C-sharp.

Could you please explain why C++ is so popular for image processing?

Thanks!

At the "dawn" of modern image processing/computer vision, C++ had the best combination of speed and industry support (meaning tools, libraries, etc and infrastructure to support them).

OpenCV is natively written in C++, look at the OpenCV about page, and you can quickly get a sense for these historical currents.


Plot Magnitude from RBJ Biquad Currently I'm working out plotting biquads in a VST/AU plugin and have problems with the graphics part. I posted in several forums but haven't got a really useful answer for this case.

First, tell me if I'm right how to get the magnitude response. For this I use this equation
http://rs-met.com/documents/dsp/BasicDigitalFilters.pdf

Am I right that I have to get an array of values, the "points" which are my Y-axis values in my filter plot?

So far, so good. Now, how should I move on to plot this magnitude on screen? I'm really on a point I don't know how to go on.

I'm using the WDL/IPlug framework in which I can use DrawPoints(x,y) and DrawLine (x,y,x1,y1) classes for drawing.

Any help would be appreciated, especially code examples/snippets (I'm learning easier when I can see it...)!

Equation (18) in the pdf is correct to calculate the magnitude response for biquad filters.

As you might know the magnitude response is plotted over the frequency. So you will need to define your points on your x-axis, which resembles the frequency $f$.

In digital systems your frequency will go up to $f_s/2$, where $f_s$ is the sample rate. The $\omega$ is just $2\pi f$. So you just have to subsitute $\omega=2\pi f$ in Eq.(18) and calculate the magnitude response $|H(e^{j 2\pi f})|$ (y-axis) for desired frequency points (x-axis).

You will end up with two vectors. One is the vector of frequencies $\mathbb{f}$ that you define. And the other vector is your calculated magnitude response $|H(e^{j 2\pi f})|$.

I have not used WDL/IPlug frameworks, but I guess you can then just plot your graph using the DrawPoints(x,y) function, using the calculated value pairs. You may have to look into the documentation of the frameworks for further hints.
Depending on the resolution of the plot, you may want to increase the step size of your frequency vector, to increase performance.


FFT peak estimation: Zero-padding vs signal repetitioon I am working on obtaining the vortex shedding frequencies behind a 2D square using CFD. Time-varying data is being extracted from a specific cell behind the square. I've written a little python script which performs FFT analysis on the signal and spits out the dominant peak frequency.

Because the shedding frequencies are so low (0.1 - 0.2 Hz) I added a zero-padding routine as well as a quadratic interpolator in order to have accuracy to several digits. My background is in fluid mechanics, not signal analysis. So my question is the following:

Would it be more beneficial to simply take my input signal (of say, 200 seconds of data) and repeat it several times (therefore increasing the length of time and subsequently the frequency resolution) or to keep the zero-padding technique?

Would the former method be much more computationally intensive? FYI, the signal is nearly sinusoidal.

Thanks

If you repeat your data, you will have discontinuities where the sequence ends meet.  These will cause spurious peaks in your spectrum.  You could apply a window to the data to taper it at the edges, but that would be a fair amount of work, and you'd need to work out a window that doesn't color your spectrum.

I would just go with zero padding.  It will mess with the amplitude, but that should be fairly easy to deal with (scale by constant factor after FFT.)


Characterization of transfer functions with no local peaks Assume that you are giving an arbitrary amplitude frequency response $A(\omega)=|H(j\omega)|$

Is there a characterization that ensures that $A$ is monotone? i.e, $A$ has a global maximum at $\omega=0$ and does not have a local peak at some other frequency?

thanks

As far as I know there is no simple characterization of filters with a monotonic magnitude response. I could only come up with a weak time domain characterization of LTI systems with monotonically decreasing amplitude for positive frequencies (cf. Eq. $(6)$ below).

Let $h(t)$ be the real-valued impulse response of an LTI system, and let $H(\omega)=A(\omega)e^{j\phi(\omega)}$ be its Fourier transform. Furthermore, define $A(\omega)$ and $\phi(\omega)$ such that $A(\omega)\ge 0$ is satisfied. Now we require

$$\frac{dA(\omega)}{d\omega}\le 0,\qquad\omega>0\tag{1}$$

Let $r_h(t)$ be the deterministic auto-correlation function of $h(t)$:

$$r_h(t)=\int_{-\infty}^{\infty}h(\tau)h(\tau+t)d\tau\tag{2}$$

Note that the Fourier transform of $r_h(t)$ is given by $A^2(\omega)$. By a basic property of the Fourier transform we have

$$\begin{align}t\,r_h(t)&=\mathcal{F}^{-1}\left\{j\frac{dA^2(\omega)}{d\omega}\right\}\\&=\frac{j}{2\pi}\int_{-\infty}^{\infty}\frac{dA^2(\omega)}{d\omega}e^{j\omega t}d\omega\\&=-\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{dA^2(\omega)}{d\omega}\sin(\omega t)d\omega\\&=-\frac{1}{\pi}\int_{0}^{\infty}\frac{dA^2(\omega)}{d\omega}\sin(\omega t)d\omega\tag{3}\end{align}$$

Note that since $A(\omega)\ge 0$, condition $(1)$ implies

$$\frac{dA^2(\omega)}{d\omega}\le 0,\qquad\omega>0\tag{4}$$

From $(3)$ and $(4)$ we get the following inequality:

$$|t\, r_h(t)|\le -\frac{1}{\pi}\int_{0}^{\infty}\frac{dA^2(\omega)}{d\omega}d\omega=\frac{A^2(0)-A^2(\infty)}{\pi}\tag{5}$$

which finally results in

$$|r_h(t)|\le \frac{A^2(0)-A^2(\infty)}{\pi |t|},\qquad t\neq 0\tag{6}$$

In sum, condition $(1)$ implies that the deterministic auto-correlation of the corresponding impulse response must satisfy inequality $(6)$.


Application of Parsevals theorem with Welch method For my project I had to manually code the welch method using the code below. Pretty much it involves finding the spectral density through fft and incorporating windowing and segmenting with overlap.

I need to know if my PSD is correct so I am checking through 2 methods.

root of area under the PSD should be the RMS (E[x^2]) and since I have taken away my mean from the data, this should equal the standard deviation of the data
Checking if Parseval's theorem holds.
There are two main issues, I am having trouble with. Firstly I am getting a 0Hz peak even after removing the mean at the start. The only way I can remove the 0Hz peak is to take the mean away after I have segmented and multiplied the window. (#### Detrend line ####)

With the peak, point 1 holds up well but without the peak the standard deviation is slightly different from the RMS but thats not a huge concern.

My second issue is that I cant get the Parseval's theorem to hold. Is there something I'm doing wrong? Is it necessary for the Parseval's theory to hold? Please let me know if my theory is wrong.

Couple of notes about the code:

apologies if it isnt efficient coding, I'm still a novice
dividing the fft by (fs*(win*win).sum()) is something I found through research and also a part of the inbuilt Welch code
lastly I am using trapezoidal rule for integration
my FFT energy is a couple of magnitudes smaller (ie E-09)

With 0Hz peak removed

With 0Hz peak

import numpy as np
import matplotlib.pyplot as plt
import os
import math
from scipy.fftpack import fft, ifft, fftfreq, rfft
from scipy.signal import detrend, get_window

plt.clf()
plt.close()

###############################################################################

config = '1perc_damping'

dir0 = os.getcwd()
dir1 = dir0+'/'+config
dir2 = dir1+'/Summary files'

data = np.loadtxt(dir1 + '/000deg_500rpm_An_Mxx.dat')

###############################################################################

fs = 625
nperseg = 1024
noverlap = 512
window_meth = 'hanning'

###############################################################################

#detrend
data = data - np.mean(data)

#get windowing method
win = get_window(window_meth, nperseg)

# number of segments
nseg = np.ceil((len(data)-nperseg)/(nperseg-noverlap)) + 1

total_n = nperseg + (nseg-1)*(nperseg-noverlap)

# padding with zeros
if total_n > len(data):
    n = total_n - len(data)
    padding = np.zeros(np.int(n))
    data = np.concatenate((data,padding))

# fft

fft_data = 0

for i in range(0,np.int(nseg)):

    start = i*(nperseg-noverlap)

    f_data = np.multiply(data[start:start+nperseg],win)
    f_data = f_data- np.mean(f_data)   #### Detrend line ###

    transf = fft(f_data)
    transf_c = (transf*np.conj(transf))/(fs*(win*win).sum())

    fft_data = fft_data + transf_c

Pxx = np.real(fft_data)/(nseg)
fx = fftfreq(nperseg, 1/fs)

plt.figure(1)
plt.plot(fx[0:512],2*Pxx[0:512])
#plt.plot(fx,Pxx)

plt.figure(2)
x = np.linspace(0,len(data),len(data))
plt.plot(x,data)

#Finding RMS
l = np.int(len(fx)/2 -1)
c = 2
A = 0

for i in range(1,l+1):
    A = A + 0.5*(c*Pxx[i]+c*Pxx[i-1])*(fx[i]-fx[i-1])

RMS = np.sqrt(A)

# Checking Parseval's theorem
data_energy = 0
FFT_energy = 0

for i in range(1,l+1):
    FFT_energy = FFT_energy + 0.5*(np.power(c*Pxx[i],2)+np.power(c*Pxx[i-1],2))*(fx[i]-fx[i-1])

for i in range(1,len(data)):
    data_energy = data_energy + 0.5*(np.power(abs(data[][1][i]),2)+np.power(abs(data[i-1]),2))*(i-(i-1))*(1/625)

RMS value

  Firstly I am getting a 0Hz peak even after removing the mean at the start. The only way I can remove the 0Hz peak is to take the mean away after I have segmented and multiplied the window.

The mean doesn't have to remain at 0 when you apply a window. Changing the weight of a particular sample would shift the mean accordingly (e.g. a weight greater than 1 on a positive sample would increase the mean). If you must remove the 0Hz peak, the only guaranteed way to do so is right before the FFT (or right after the FFT, but that second option wouldn't remove the spectral leakage to adjacent bins).

Parseval's relationship

  My second issue is that I cant get the Parseval's theorem to hold. Is there something I'm doing wrong? Is it necessary for the Parseval's theory to hold? Please let me know if my theory is wrong

Since you are dealing with discrete time signals, you should be using the specific formulation of Parseval's theorem for discrete time signals. This formulation is derived directly from the definition of the Discrete Fourier Transform. It then has to hold between the input and output data to the fft call. Note however that different implementations may use different normalization which can affect the exact formulation by a constant scaling factor. With the normalization used in numpy.fft.fft, you should have the following relationship:

$$
  \sum_{n=0}^{N-1} |x_n|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |X_k|^2
$$

Note however that your code currently incorrectly throws away the imaginary parts of your computed spectrum on the following line:

Pxx = np.real(fft_data)/(nseg)

You final computation should instead make use of squared-magnitudes. Since you are later taking the square of each frequency bin value, you may take the absolute value:

Pxx = np.abs(fft_data)/(nseg)

Your will also need to remove the frequency step part in your final integration of FFT_energy and stick with the straight summation given in the discrete time formulation of Parseval's theorem given above.

for i in range(0,len(Pxx)):
    FFT_energy = FFT_energy + np.power(Pxx[i],2)
FFT_energy = FFT_energy / len(Pxx)

for i in range(0,len(data)):
    data_energy = data_energy + np.power(abs(data[i]),2)


Extract hand written text from US census images As many will know, images of all US censuses going back to 1790 are available online http://blog.archive.org/2010/12/22/1790-1930-u-s-census-records-available-free/. Unfortunately the quality of the images is very poor in places and, as far as I know, there are no free processed versions of them available.  If you take an example like the following from the 1920 census, what would be a good method to extract relatively clean images of the handwritten text from the background image and noise?

The original pdf for the page is at
http://wikisend.com/download/737634/page.pdf  as I don't have the reputation points to post it here yet (can someone with more rep upload it?).

I couldn't open your pdf (antivirus said it was suspicious) so i downloaded an image myself, I don't know if its the right one, but its also from a 1920 census

You could perform image dilation on the document until the text is gone. You can do multiple passes if needed. I just did one. We know dilation will make shrink the black sections (or the text) this essentially gives us just the background.

We then remove the background from the image. This can be tricky depending on the language (I used matlab, and had to make sure to keep the right range and units for the image to display properly) then you have a cleaner image. You should be able to run a ocr engine on the result

and here is a zoomed in section just so you can see the text quality doesn't degrade much with this operation

and of course the code. As I said it is done in matlab but openCv(in many languages), python, or just about any other image library will be able to perform a similar operation

im = rgb2gray(color_im);
bg = imdilate(im,strel('disk',3));
%we do a complement becasue the bg subtraction makes the text white and
%the background black. so we invert it
clean = imcomplement(abs(bg-im));

subplot(1,3,1); imshow(im); title('original')
subplot(1,3,2); imshow(bg); title('after imdilate with 3 pixel circle')
subplot(1,3,3); imshow(clean); title('original minus background')

and here is the image I used. Its the first page of the 1920 census https://archive.org/details/fourteenthcensus02unit

edit

I forgot you said it was handwritten images. so I got another sample, same exact code, new image

notice the detail even in very dark obscured places, as the zoomed in view shows. You would need to run some other filters, but removing the main background noise is accomplished using dilation


Meaning of sampling frequency? Say there is a communication system, which can manage sampling frequency up to 4 MHz (given by software). Can the system be able to process signals with maximum frequency more than 2 MHz? Because according to sampling theorem sampling frequency > 2*maximum frequency.

Probably no, but... maybe. The classical sampling theorem is often misstated: if a continuous signal is band-limited, it can be recovered exactly (in theory) if the sampling frequency is twice its maximum frequency.
However, there are cases where signals can be recovered above the Nyquist limit, provided that we have more information, because ambiguities may happen  in the recovery:

finite-time signal cannot be strictly band-limited, but processing can afford a small amount of distorsion,
some band-limited signals can be recovered with about twice their bandwidth (and not just twice the max frequency).
there are more generic sampling theorems (like: do you sample the derivative of the signal at the same time), but more involved as well.

Generally, the above possibilities are not implemented in traditional hardware or software.


How to determine the step response given a transfer function? I'm new to Signal Processing. I would like to know how to draw the step response given a transfer function.

For example, given $G(s)=\frac{(1-\frac{s}{3})}{{(\frac{s}{0.1}+1)(\frac{s}{100}+1)^2}}\cdot 10$ a transfer function, I calculated the initial value with the "Initial Value Theorem": $g(0^+)=\lim_{s\rightarrow \infty}sG(s)=0$, but other than that, I would also like to calculate the final value and the settling time.

I've been also given the Bode diagrams (I don't know if it helps):

I think $G(s)$ is a transfer function of the third order, so do the second-order laws also apply there, where the settling time is calculated as $t_{s}=\frac{4}{\zeta \omega_{n}}$?

The short answer is that the step response is the inverse Laplace Transform of $\frac{1}{s}G(s)$.

Here are some hints that should help in being able to solve the problem. Since it is likely a homework problem I wanted to provide more details as relevant background and steps toward a solution rather than detail the mathematics for the specific solution. A plot of the resulting step response is included at the end to validate the solution.

First, a review of the simple case of real negative poles on the s-plane along with their inverse (unilateral) Laplace Transform.

As shown, a pole given by the transfer function $H(s) = \frac{1}{s+\alpha}$ has an inverse Laplace transform $x(t) = e^{-\alpha t}$ for $t\ge 0$. This is easy to prove but is also well tabulated in all tables of the Laplace Transforms such as at the bottom of this page: https://en.wikipedia.org/wiki/Laplace_transform

Specifically if $H(s)$ describes the transfer function of the system (output divided by input when the units are $s$), $h(t)$ is the impulse response of the system. As an aside demonstrating one utility of the Laplace Transform, if we restrict $s$ to just be the imaginary axis, we get the frequency response $H(j\omega)$ which in this case would be a low pass filter. Whether $s$ is real or imaginary, we typically refer to the s-domain as the "frequency domain" and all points on the surface are given as pole locations by the transforms of $e^{st}$. So with $h(t)$ being the impulse response, what we see as graphed is the response in time at the output if an impulse was presented at the input. Notice that the pole closest to the imaginary axis takes the longest to decay in its impulse response. This is a key take-away and is the reason we call the negative poles closest to the axis "dominant poles". Like an unwanted guest, they take way too long to go away. Poles further away from the imaginary axis decay much more rapidly in time.

Now consider the summation of these two transfer functions:

$$G(s) = H_1(s) + H_2(s)$$

A summation of two systems in the frequency domain is a parallel process in their signal flow implementation as depicted in the figure below. In contrast the function $H_1(s)H_2(s)$ would be implemented as the cascade of the two systems in series.

Now for the system above represented by the summation of our two systems, the total output due to an impulse at the input would be the sum of the impulse response from each system. Here we see how the pole that was further away from the imaginary axis quickly decays so does not significantly contribute to the final settling time of the output. This settling time (and the overall shape of the response if the poles are far enough apart) is dominated by the pole closest to the imaginary axis, our dominant pole.

Further this combined system is given as:

$$G(s) = \frac{1}{s+b} + \frac{1}{s+a} = \frac{2s+a+b}{(s+a)(s+b)}$$

Notice how adding them together introduced a zero (at $z_1=-(a+b)/2$) but did not change the pole locations. Importantly this shows how given the general combined expression on the right side of the above equation, the process of Partial Fraction Expansions can be used to put this in the form of the parallel (summed) system as described by the left side of the above equation. This is then very easy to solve for the inverse transform,  providing the impulse response, here given as:

$$g(t) = e^{-bt} + e^{-at}$$

Step Response

The above shows how to solve for the impulse response, which is the inverse Laplace Transform of the transfer function $G(s)$.  To get the step response note that the Laplace Transform of an integrator is $1/s$ (also in the Wikipedia link above and a very important takeaway). Similarly the Laplace Transform of a differentiator is $s$ (less the initial condition). Yes commit these to memory as they will come up again and again and again. So imagine if we prefaced our system under test with an integrator: The integratiion of an impulse is a step! So for any given system, if we simply multiply it's transfer function by $1/s$ (which means putting an integrator in cascade or series with the system), the output defined by the inverse Laplace Transform of that result will be the step response! It's that simple. Taking that further if we multiplied by $1/s^2$ we would get a ramp response, etc.

Solving this again is aided by Partial Fraction Expansion:

$$\frac{1}{s}H_1(s) = \frac{1}{s(s+b)}= \frac{1}{b}\bigg(\frac{1}{s} -\frac{1}{s+b} \bigg)$$

From which we can easily derive the inverse Laplace Transform providing the step response $s(t)$:

$$s(t) = \frac{1}{b}\mathscr{L}^{-1}\bigg(\frac{1}{s} -\frac{1}{s+b} \bigg) = \frac{1}{b}(1-e^{-bt})$$

To easily see how the consideration of the dominant pole still applies, consider that a step at the input to the system is an infinite series of impulses each at time $t$ (and to be accurate each is weighted by $dt$ but that is just a scaling). Each impulse produces it's related impulse response, with the responses from llater impulses (accumulating) with the tails of the responses from previous impulses. The accumulation of the weighted impulse responses is described mathematically as convolution but results in the step response at the output. Therefore the dominant pole would also dominate the resulting step response.

With that background in mind, here are the main points specific to solving the OP's problem:

The step response can reasonably estimated from the dominant pole alone.

For an accurate and complete solution, use Partial Fraction Expansion to factor decompose $G(s)$ into three parallel systems $G_1(s)$, $G_2(s)$, $G_3(s)$, which can then each be easily solved for their inverse Laplace Transforms (this would provide the impulse response of $G(s)$.

The step response for each can be solved as described above by first multiplying by $1/s$ and then solving for the inverse Laplace Transform. The result for this is plotted below with $G_1(s)$ containing the dominant pole $s=-0.1$, $G_2(s)$ containing the weaker pole  at $s=-100$ and $G_3(s)$ containing the two repeated poles at $s=-100$ (consistent with Partial Fraction Expansion for repeated poles).

It can be confirmed how closely $G_1(s)$ follows the step response for the system with only the dominant pole given by $1/(s+0.1)$. Also notice the significant difference in time and amplitude scales between the three plots. $G_1(s)$ takes over 50 seconds to settle while $G_2(s)$ and $G_3(s)$ are nearly settled at 50 milliseconds. The sum of all three of these will closely match $G_1(s)$ alone. Further the final values for $G_2(s)$ and $G_3(s)$ are significantly smaller than the final value for $G_1(s)$, but the sum of all three final values will match the final value for combined system $G(s)$.

Comparing the overall step response to that given by the dominant pole highlights the significance of recognizing the dominant pole and how that can be used to simplify most practical problems (I have time now to write this all out from all the time I saved before by recognizing the dominant pole!):

The more significant difference between the two is in the very beginning, so zooming in to show that detail:

The MATLAB/Octave with the control toolbox to plot the step response for $G(s)$ is:

num = -[1 -3]*10000/3;
den = conv([1 0.1], [1 200 10000]);
G = tf(num, den)
step(G);


Why does resampling change the frequency range? I am currently working on a set of .wav files which, unfortunately, do not have the same sampling rate - I need to resample them for the purpose of feeding them to a neural net.
So - using Sox I resampled every file - with the command:

sox -v 0.99 filein -G -r 192000 fileout

However, having a look at the resampled files: ,
the frequency range has changed. Which is confusing me. I have put several hours into finding out why this is happening but i can't seem to find the reason.
Can somebody explain?

EDIT1: Just realized that this example is misleading as I am referring to the change in maximum Frequency, (the y-labels)

Uniformly sampled pulse code modulation (PCM) audio such as WAV files only encodes frequencies up to the Nyquist frequency, which is half the sampling frequency. So for a sampling frequency of 500 kHz you can have frequencies up to 250 kHz, and for a sampling frequency of 192 kHz you can have frequencies up to 96 kHz. This is reflected in the range of the Y labels.


Subsampling property of wavelet transform One of the properties I have seen for isotropic tight wavelet frames is

$$\sum_{i\in\mathbb{Z}} \left|h(2^i\omega)\right|^2 = 1$$

where $h(\omega)$ is the frequency spectrum of the original wavelet. See page 8 in A unifying parametric framework for steerable wavelet transforms.

I have always assumed that this enables the downsampling often seen in wavelet decompositions, which start off with a wavelet satisfying $h(\omega) = 1\, \forall \omega > \pi/2$, meaning that the top half of the spectrum is covered and thus subsampling of the low pass remainder is possible because there are no longer any frequencies above $\pi/2$.

Is this correct?
Or can you perform downsampling still having frequencies above $\pi/2$?
If one does not need downsampling or wants to have different partitions, is the following condition sufficient,
$$\sum_{i=1}^{I} \left|h_i(\omega)\right|^2 = 1$$
where $\{h_i(\omega)\}$ are the frequency responses of the wavelets in the frame?

I am not completely confident with the full theory of wavelet frames, however I will share some pointers and my beliefs, do not take any of  the following for a truth. The property you show is (partly) a characterization of an orthogonal system. In Hernandez & Weiss, A first course on wavelets, 1996, you find Theorem 2.12:

If $\psi\in L^2(\mathbb{R})$ is a band-limited function such that $\hat{\psi}$ is zero in a neighborhood of the origin and $\{ \psi_{j,k}: j,k\in \mathbb{Z}\}$ is an orthonormal system, then this system is coomplete if and only if
$$\sum_{j\in \mathbb{Z}} |\hat{\psi}(2^j \xi)|^2=1 \quad \text{a.e. on }\mathbb{R}-\{0\}$$
and
$$\sum_{j=0}^{\infty} |\hat{\psi}(2^j \xi)\hat{\psi}(2^j (\xi+2k\pi))| =0\quad \text{a.e. on }\mathbb{R},\quad k\in 2\mathbb{Z}+1$$
with the system defined as:
$$ \psi_{j,k}(x)= 2^{\frac{j}{2}}\psi(2^j x-k)\,.$$
So the conditions are a little bit stricter for an orthonormal system.
So it seems rather to me a characterization in 1D (f band-limited wavelets), than can help define radial 2D wavelets, than a characterization of radial wavelet per se. So if I understand your first question well, this imply the possibility of discrete wavelets and downsampling.

Since it is a characterization, I (very unsure) would say that you can downsample, but you have no guarantee about the completeness. I do not understand your condition  $$h(\omega)=1, \forall \omega>\pi/2, h(\omega)=1, \forall\omega>\pi/2,$$
but I can add that in Chapter 3.4 of the above book, you have characterization of wavelets with support in $[-\frac{8}{3}\pi,-\frac{2}{3}\pi] \cup [\frac{2}{3}\pi,\frac{8}{3}\pi] $, so I would say yes, as $\frac{8}{3}> \frac{1}{2}$. I would say that the limit on $\pi/2$ is not an issue because of the dilations of the mother wavelet.

For the last question, without further information on the $h_i$, it is unlikely that you generate a wavelet system (possibly because of refinement equations), but as long as you cover the whole spectrum, who really cares? You get a decomposition that is energy preserving, but I am not sure yet it suffices to define a tight frames.

Additional lectures:

Affine Density, Frame Bounds, and the Admissibility Condition for Wavelet Frames, 2007, Gitta Kutyniok, with results in frames with generalized samplings  $ \psi_{j,k}(x)= a^{\frac{j}{2}}\psi(a^j x-b)$
A first course on wavelets, 1996, Eugenio Hernandez and Guido L. Weiss
Ten Lectures on Wavelets, 1992, I. Daubechies, esp. Chapter 3
Tight Wavelet Frame Construction and its application for Image Processing, 2005, Kyunglim Nam
Frames and bases. An introductory course, 2008, Christensen, O.


I tried two approaches and gained the different conclusions of judging the stability of the transfer function of the system We want to judge whether the system is stable or not.
Given the below transfer function.
$$ H\left( z  \right) =\frac{\left( 1+2 z^{-1}  \right) }{\left( 2+z^{-1}  \right) }   $$
$$ H\left( z  \right) =2-\frac{3}{2+z^{-1} }   $$
$$  2+\frac{-3}{2+z^{-1} } $$
$$S:=\frac{-3}{2+z^{-1} }$$
$$ S ~~\text{is the sum of each element of the geometric sequence.}  $$
$$ -3 ~~ \leftarrow~~ \text{initial term}  $$
$$ z^{-1} ~~ \leftarrow~~ \text{common ratio}  $$
$$ i \geq1 \rightarrow \text{ith term} = -3 \cdot \left( z^{-1}  \right) ^{i-1}  $$
$$ = -3 \cdot z^{1-i} =-3 \cdot z^{-0}~,~-3z^{-1}~,~-3 z^{-2} ~,~ \cdot\cdot\cdot     $$
$$ H\left( z^{}   \right) =\sum_{ n=-\infty  }^{ \infty  } h\left[ n \right]  z^{-n}   $$
$$ \sum_{ n=-\infty  }^{ \infty  } \left( 2 \delta\left[n  \right]  +\left( -3 \right) u\left[ n  \right]  \right)z^{-n}     $$
$$\displaystyle \therefore ~~ h\left[ n \right] =2 \delta\left[n  \right]  +\left( -3 \right) u\left[ n  \right]  $$
$$  \sum_{ n=-\infty  }^{ \infty  } \left| h\left[ n \right]  \right|    $$
$$ = \sum_{ n=-\infty  }^{ \infty  } \left|  2 \delta\left[n  \right]  +\left( -3 \right) u\left[ n  \right]\right|  = \infty  $$
$$  \therefore ~~ ~~\text{The system is unstable.}~~   $$
However from the another approach,
$$ H\left( z  \right) =2-\frac{3}{2+z^{-1} }   $$
$$ = \frac{2 \left( 2+ z^{-1}  \right) -3}{2+z^{-1} }   $$
$$2+z^{-1}  = 0  $$
$$  2+\frac{1}{z^{} } =0  $$
$$  \frac{1}{z^{} } =-2  $$
$$  z=-\frac{1}{2} =-0.5 ~~ \leftarrow~~ ~~\text{pole}~~   $$
$$  ~~\text{Since }~~ \left| \text{pole}  \right|  =\left| -0.5 \right| =\left| 0.5 \right| <1 ~~\text{is held, the system is stable.}~~  $$
Why the different concolusions were gained?
What I've been missing?

The system is stable. Your mistake is in the application of the formula for the geometric series:
$$\begin{align}-\frac{3}{2+z^{-1}}&=-\frac32\frac{1}{1+\frac12z^{-1}}\\&=-\frac32\sum_{n=0}^{\infty}\left(-\frac12\right)^nz^{-n}\end{align}$$
Hence, the system's impulse response is
$$h[n]=2\delta[n]-\frac32\left(-\frac12\right)^nu[n]$$


Digital filter design of time series for specified frequency response function I am currently working with vibration measurements in structures. In the netherlands there is a guideline for verifying vibration measurements for damage to machinery. This is the so-called "SBR Trillingsrichtlijn". In this guideline a frequency weighting function is specified to modify the time series from vibration measurements:

with small f being the frequency in Hz, f_0 a constant frequency of 5.6 Hz and v_0 a predefined constant velocity of 1 mm/s. To illustrate the filter behaviour i show it graphically here as well (axes are linear):

In post-processing this is fine as i can perform an FFT on the time series data, apply this weighing function in the frequency domain (through multiplication with the FFT transformed time series) and then do an inverse FFT to get the adjusted time signal data.

My question is the following: the equipment we use to measure vibrations does this weighing internally using a digital filtering on the time series in real-time. To verify the data i would like to design a digital filter corresponding to the weighing function above, can anyone point me to a resource how to approach this? From what i have looked up so far on IIR / FIR filter design usaully the examples start from a transfer function in the s-domain or z-domain (which correspond to a laplace transform in discrete or continuous time) but now i only have a frequency response in the fourier / frequency domain.

EDIT: Based on suggestion of @Hilmar it appears a simple first order high pass function. i made a python implementation to check the method of FFT and the highpass filter. shown in the graph here. Strangely the frequency response of the highpass filter and SBR specified response match exactly but the implementation using the highpass has some artefacts whereby non existant frequencies are generated with non-zero amplitudes:

this can be seen in the green dotted line which is a mock signal filtered in python using the scipy package.

This is the classical problem of filter design:

  You've got a frequency response that you want, how to implement it using a FIR?

I'm not going to lay out all the theory here, because it can easily be found by looking for Fourier approximation methods for FIR design or windowing methods, but the idea is easy enough:

Take your frequency response. Sample it at many points.
Take the inverse FT. You get something long, which, when convolved with a signal, would have the sampled frequency response (because of the convolution theory of the Fourier Transform).
Because the result is a very long filter and exhibits undesirable fringe behaviour, apply a window function to cut off some of the coefficients. The result is your filter taps.

There's iterative methods of increasing the match between desired and actual response given a limited amount of freedom by a limited number of taps; but this would exceed the scope of the question.


Impulse response of system with $\max\left\{\lvert x[n] \rvert\right\}\ge\max\left\{\lvert y[n]\rvert \right\}$ Find a relation/condition for the impulse response, so for any $x[n]$ (input) with $y[n]$ (output) this relation is true:
$$ \max\left\{\lvert x[n] \rvert\right\} \ge \max\left\{\lvert y[n]\rvert \right\}$$
The system is discrete, linear and time invariant.

I tried to replace the $y[n]$ with the convolution formula,but it doesn't work. I tried to use Fourier transformation, but it doesn't work either. Any suggestions?

Not sure what exactly the question is: find some LTI system that meets this constraint (which is trivial) or create a testable condition for any LTI system. (Which is not trivial).
For the latter: Discrete LTI systems are guaranteed to not increase the maximum amplitude if the absolute sum of the impulse response is smaller than $1$, i.e.
$$\sum |h(n)| \leq 1 $$
The proof is a bit tedious, so I'll skip it unless there is need for it.


Difference between filters in serial and parallel I am trying to understand the difference between filters in serial and filters in parallel. As I understand it, serial is the more natural approach: A is applied then B, and the result is the product. Parallel filters confuse me though. If I think about it, I'd assume their outputs are added. Is this like a resistor in an electric circuit?
The book I'm reading (DSP by Jonathan Stein) says capicitors in series are high pass, while in parallel are low pass. Inductors the other way around. Not knowing much about electrical engineering, this confuses me.
Is the understanding I have shared at least correct? Any help with the capacitor and inductor example? Thanks!

If you're talking block diagrams, a pair of filters in series would yield
$$H(z) = H_1(z) H_2(z)$$
The only sensible way to connect filters in "parallel" in a DSP application would be to sum them:
$$H(z) = H_1(z) + H_2(z)$$
The response would be completely different (work it out with a couple of different transfer functions for $H_1$ and $H_2$ and see what I mean).
Analogies should be an aid to understanding, not a barrier.  That book must be assuming folks who are familiar with circuits.  So I'm going to give you a short course -- but finding a book that doesn't assume an intimate knowledge of circuits may be a good idea for you.
Here's two lowpass filters and two highpass filters, and two filters usin resistors and caps, and two filters using resistors and inductors.
Basically, a capacitor is an open-circuit at low frequencies, and a short at high frequencies.  So when you arrange the capacitor as a shunt to ground, at high frequencies all the signal goes to ground, but at low frequencies it all goes to $V_{out}$.  Conversely, when you use a cap as a series element and a resistor to ground, at low frequency the cap blocks the signal, and the resistor passes whatever is left to ground.
An inductor is the opposite -- it's a short circuit at low frequencies, and an open circuit at high frequencies.  So the circuits are opposite.
But I want to stress here -- if what I just said still confuses you, get a different book.  DSP is about math, not electrical circuits.  If the author is relying on an intuition about electrical circuits to give you an intuition about DSP, and you don't have an intuition about electrical circuits, then you need a book that's done a better job aiming itself at you, not the electrical engineer standing next to you.
Even if it's for a class, you may want to get a different book as a suppliment.


What Is the Significance of a Large Residual When Applying Principal Component Analysis? I am using Matlab function PCA (principal component analysis) to reduce the dimensionality of a data set with approximately 20 000 observations x 100 dimensions.

After having obtained the principal component coefficients of the data I recreated the input signal in the original coordinate system using the transformation matrix from the PCA function. This yielded a very large residual when comparing with the input signal. I have tried around with different data sets and sizes and it appears to be commonplace to have a large residual. I am not sure yet whether it is due to round-off errors or high SNR in the input data. The dimensionality reduction could of course still be useful, but is this something that one should be cautious about when performing principal component analysis? Or is there another metric that is better to assess the performance of the principal component analysis?

PCA dimensionality reduction seeks for a linear transformation, mapping the data to a lower dimensional space. The linearity comes from the mapping, which is simply a rotation (and translation if you take into account the de-meaning). In the end, it is nothing but a projection of the data.

So if your low dimensional embedding is not reconstructing the data truthfully, then this means that your data is not suitable for such linear mapping, i.e. not linearly separable. You should look into other techniques such as linear discriminant analysis (canonical variates), Kernel-PCA, which applies kernel trick to achieve non-linearity or manifold learning.


log-likelihhood function for N sample of data if $x(t)=b A e^{ j\omega t} + e(t)$ for $t= 1,2,...,N$ where $b$ is a parameter, $A$ is a vector $M \times 1$, $e(t)$ is a white Gaussian noise with covariance matrix of $Q$ theh what is log-likelihood function for that?

I’m assuming that the trace is what is confusing, so let me give you some clues.

First an identity $$ Trace(\mathbf{ABC)})= Trace(\mathbf{BCA}) $$

and for scaler $x$,
$$ Trace(x) =x$$

Let $\mathbf{a}$ and $\mathbf{b}$ be $N \times 1$ vectors and $\mathbf{C} $is a $N \times N$ matrix, so $$ \mathbf{ a^T C b} = Trace ( \mathbf
{C \,b a^T}) $$

The rest you should be able to figure out using the definition of a Gaussian, independence, and properties of the log.


Deep zeros in the spectrum of the input data? I wanted to know the meaning of "deep zeros" in the spectrum of the input data. I came across this terminology while going through the introduction of chapter 2 (Bussgang techniques for blind deconvolution and equalization) of the book named "Blind Deconvolution" by Simon Haykin.

Thanks,
JK

"Deep zeros" in the input spectrum are caused by a channel which strongly attenuates one or several frequencies, often due to multipath fading. This causes problems for a (linear) equalizer because it tries to compensate for these "deep zeros" by strong amplification of the corresponding frequency bands, which leads to noise enhancement. Non-linear equalizers are better at dealing with such problematic channel conditions.


As far as I can see, there are two little mistakes. First of all, you should use dct2() instead of dct(). And second, you shouldn't use the 'transpose' function because you want the DCT matrix not only transposed but also conjugated. So you should do the following:
$$\tt{T=dctmtx(4); DCTB = T*B*T'}$$
And this should equal $\tt dct2(B)$.

EDIT: As correctly pointed out by Peter K., the DCT matrix is real-valued, so indeed the 'transpose' operation and the ' operation give the same result. So I only noticed ONE little mistake ...


Precompute parts of an FFT on streamed data inputs I've implemented Radix-2 & Radix-4, DIT & DIF variants of the cooley-tukey FFT algorithm and they perform well on test data.

Now for a project, my algorithms will be applied to a real-time data stream of complex values. My current approach is window-based, meaning that I'm creating a buffer that stores values in a 2d array for a certain amount of time and after that, the FFT is applied to the matrix stored in the buffer.

As far as my understanding of the algorithm goes (having implemented it and read some papers on DFT/FFT in general but no prior experience in DSP), there is no way of speeding this up by doing some computation live as the data points arrive in the stream connector?

I was wondering if this is correct since I did not find any literature/paper that does this. Any hint is appreciated.

Edit: Example for clarification

Say my input is a Matrix of size N x N of complex numbers and I'm getting the values one point at a time from the data stream, can I do anything more than just buffering them in a 2d array and doing the FFT after a certain time.

A sliding DFT can be used to compute the equivalent of a new FFT as each input element arrives.  However, the cost of a sliding DFT is O(N^2), not O(NLogN).


How is the lowpass to bandpass transformation derived? I'm using a lowpass to bandpass filter transformation:

$$ s \leftarrow Q \left( \frac{s}{\omega_0} + \frac{\omega_0}{s} \right) $$

Based on the wiki article: http://en.wikipedia.org/wiki/Prototype_filter which basically is a change of variable. But I was wondering how this transformation is derived?

The reason being I need to create my own transformation so I can combine a low pass with a shelf using a single transformation from a low pass (to use in the Nth order butterworth equation). Knowing how the lowpass to bandpass transformation is derived would take me a step further to experimenting combining other filters.

In the $s$-domain, the LPF-to-BPF transformation doubles the order of the filter.  that is because the LPF has one transition from passband to stopband, but the BPF has two such transitions.

Remember that when we do frequency response, we substitute $s=j\omega$.  That is, we evaluate the $s$-plane transfer function $H(s)$ as $H(j\omega)$.  The magnitude of $H(j\omega)$ or $|H(j\omega)|$ is the gain of the filter when the input frequencey is $\omega$.  Also keep in mind that the reciprocal of $j$ is its negative: $ \frac{1}{j} = - j $.

In the LPF, the middle of the passband is at DC, $\omega=0$.  The transformation to BPF is such that maps, using a nice continuous function, the middle of the BPF passband at $\omega=\omega_0$ to $\omega=0$.  A simple, mathematically continuous mapping as such is:

$$ \omega \leftarrow A ( \omega^2 - \omega_0^2 )  \quad A \ne 0 $$

So whenever the $\omega$ on the right (which is the BPF $\omega$) is equal to either $+\omega_0$ or $-\omega_0$, the $\omega$ on the left (which is the LPF $\omega$) is zero.

now suppose we set $A \triangleq \frac{Q}{\omega_0 \omega} $ which is still non-zero (and non-infinite for $\omega \approx \omega_0$).  now you have

$$ \omega \leftarrow Q \left( \frac{\omega}{\omega_0} - \frac{\omega_0}{\omega} \right) $$

$$ \begin{align}
j\omega & \leftarrow Q \left( j\frac{\omega}{\omega_0} - j\frac{\omega_0}{\omega} \right) \\
& = Q \left( \frac{j\omega}{\omega_0} + \frac{\omega_0}{j\omega} \right)  \\
\end{align} $$

or

$$ s \leftarrow Q \left( \frac{s}{\omega_0} + \frac{\omega_0}{s} \right)  $$

so this means, regarding frequency response:

$$ \begin{align}
H_\text{BPF}(j\omega) & = H_\text{LPF}\left(j Q \left( \frac{\omega}{\omega_0} - \frac{\omega_0}{\omega} \right) \right) \\
& = H_\text{LPF}\left(Q \left( \frac{j\omega}{\omega_0} + \frac{\omega_0}{j\omega} \right) \right) \\
\end{align} $$

or, regarding transfer function:

$$ H_\text{BPF}(s) = H_\text{LPF}\left(Q \left( \frac{s}{\omega_0} + \frac{\omega_0}{s} \right) \right) $$