Dominant Orientation using Fourier Power Spectrum I have been using Fourier Transform Power Spectrum for estimating the dominant direction in some images.

However, I am now working on extremely small images of size 8X8 pixel and I have found that even at this level there is some result regarding the dominant orientation of the image. As the image size is 8X8, the total number of pixel will be 64 and as shown in the figure below, there will be less number of angles.

I have simply done polar conversion  of the Fourier Transform Power Spectrum of the image which will be of same size (8X8), then I have taken the summation of the intensity for each angle and plotted them against the angle.

Now, from the plot I can find the maximum angle but I want to increase the accuracy of my result. I will first find the dominant direction of an image (a medical image) and then I will find the changes in the diseased state.

Thank you so much for your help.

Subtract from the 8x8 image the mean value of its pixels, which should not affect the dominant orientation. This will remove artifacts due to windowing in the next step.
Window the image with a circularly symmetric window function that eliminates artifacts aligned with the image axes. Make sure you have chosen a window function that is sufficiently smooth to be sampled with a 8x8 grid without frequency domain aliasing.
Pad the image to 32x32 with zeros to increase the image size and thus the discrete Fourier transform (DFT) resolution. This makes the DFT result smoother and thus easier to interpolate.
Calculate for each angle an integral of the squared magnitude of the DFT result over a line corresponding to the angle. All angles must be treated equally; use same length lines. 4x4 piece-wise bicubic polynomial interpolation should suffice if you approximate the integral by a sum of point samples at the axis-aligned resolution of the DFT.

If you want to make the computation faster you can take advantage of that the integral will change quite smoothly as function of the angle, so you can first sample at the coarse resolution of at least $\text{atan}(1/8)\times 360\unicode{xb0}/(2\pi) \approx 7\unicode{xb0}$ to find where the highest peak might be, and then sample more finely around that angle.


Carrier frequency offset and carrier leakage Recently I study frequency Offset estimation due to mismatch in carrier frequency from the oscillator at the receiver. I understand how this cause problem can cause a phase shift. But I also saw that this cause problem called CARRIER LEAKAGE. So I look online and I saw the following:

Certain transmitter implementations may cause leakage of
the center frequency component. Such carrier Leakage can
occur in some transmitter because of DC offset. This issue
is seen on the receiver side as energy in the transmit center
frequency

Does carrier leakage mean an increase in interference power?? I thought that carrier offset happens at receiver and not transmitter. I also saw that carrier leakage can cause increase in the Error Vector Magnitude (EVM), I can't see the relationship between carrier leakage and EVM.

Looking forward for your discussion.

I think you are confusing two different things.

Carrier offset is the difference between the transmitter's and the receiver's local oscillator (LO). The receiver needs to eliminate or reduce the offset if it hopes to receive data consistently. This is done by estimating the offset and adjusting the frequency of the LO.
Carrier leakage causes a DC offset in the receiver. Consider the receiver mixer: it has two inputs, one is the received signal, and the other is the local oscillator. In practice, the LO signal will be present at both inputs, because it is very difficult to isolate high-frequency signals in a PCB. This causes the LO to mix with itself, producing a DC term at the mixer's output.


Evaluate the Z Transform Evaluate the Z transform of $x[n] = n^3$ where the signal is two sided.

I have tried using the basic definition of the Z transform ie.,

$$X(z) \triangleq \sum_\limits{n=-\infty}^{+\infty} x[n] \, z^{-n} $$

But what I get obviously is an infinite unbounded sum of sequence. Does that mean Z transform does not exist? IS there any alternate methods to obtain the transform? Please help.

[My path was erroneous, I made a confusion, but kept a modified answer for one-sided $\mathcal{Z}$-transforms]

With $\mathcal{Z}$-transforms, formulae with powers of $n$ often aim at making you handle basic operations on $\mathcal{Z}$-transforms: linearity, time-product/differentiation, integration, time-shift, see for instance: Properties of Z-Transform:

$$ \mathcal{Z}(nx[n]) = -z\frac{d}{dz}X(z)\,,$$
and
$$ \mathcal{Z}(x[n-n_0]) = z^{-n_0}X(z)\,.$$

Hint: if your case were a one-sided $\mathcal{Z}$-transform, one could have focused on derivation, and on the simple geometric series:

$$\sum_{n=0}^{+\infty}z^{-n} = \frac{z}{z-1}\,.$$

As of today, as pointed out by Matt and Robert, I could not make sense of the
$$\sum_{n=-\infty}^{+\infty}z^{-n}\,,$$
and I can see but a zero radius of convergence.
For what it's worth, I will anyway try to add a common trick to try to derive something. One can use a stabilization with  some $0 < a <1$:

$$\sum_{n=-\infty}^{+\infty}a^{|n|}z^{-n}\,,$$
with the secret thought that, if $a$ is just close to one, we could get almost what we wanted.
I may have made mistakes again, so far this yields the formal form:

$$ \frac{(1-a^2)z}{(z-a)(1-az)}\,.$$

This still have a singular ROC, but one can play with  $\mathcal{Z}$-transform properties.


Extracting direction of dominant sinusoids from an image in terms of angles Here is a sample image.

What I am looking to do is compute the dominant "wrinkle" directions.
The dominant wrinkle directions can be seen in the following FFT, roughly, by drawing a line through the highest valued set of points in the spectrum.

I have been told that a potential way to find this direction would be to look at the auto-correlation scipy.signal.correlate2d.
I have not worked with this tool before, so I want to make sure my idea is correct.
The way that I think this will work is by using the FFT data to find the largest components of the signal.
Say for example, I look at the largest two components.
I can construct two images, one for each of the largest components, lets call the images dominant1 and dominant2.

For each image, I could then look at:

scipy.signal.correlate2d(image,dominant1)
scipy.signal.correlate2d(image,dominant2)

I believe that this will return a matrix that shows where the features are in the image that most agree with the signals dominant1 and dominant2.

If this interpretation is correct, I am still unsure of how to extract the dominant_angle1 and dominant_angle2 say. Suppose that the horizontal direction in the image as shown above is x and the vertical is y. I would like dominant_angle1/2 to be the angle off of the x-axis.

Any comments, suggestions for better approaches, will be appreciated.

The following should work

Perform FFT

$$F(u,v) = \mathcal{F}[f](x,y)$$

Multiply magnitude by 2nd order circular harmonic and sum

$$c = \sum_u \sum_v |F(u,v)|  e^{i 2 \phi}$$

where $\phi = \arg (u + iv)$

Angle is half of argument

$$\theta = \arg(c) / 2 $$

Here is some matlab code:

I = rgb2gray(imread('eNb2pye.jpg'));
F = fft2(I);
[~,~,r,th] = FFTMesh(size(I));
c = sum(sum(abs(F) .* exp(1i.*2.*th)));
direction = angle(c)/2

which uses the following function to general the spectrum u,v and radial coordinates

function [ux,uy,r,th] = FFTMesh(sizeM)
rows = sizeM(1);
cols = sizeM(2);
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);
% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);
end


Finding Signals when the baseline varies greatly by signal set I have milliwatt consumption data from home appliances.  I'm collecting this data in hundreds of households.  Based on the machine model, the baseline value for a given household can vary from 0 to 500 milliwatts.  The baseline non-zero values vary by 5 to 15 milliwatts and the data is reported every 13 seconds.  The vast majority of data is baseline, and there is only signal when people are running the appliance.  Perhaps 10 hours a week is real data.  The signals are very pronounced, going from baseline to values exceeding 100,000.  I have 1/3 of a billion records (growing by several million records a day) so signal detection must be automated.

I really just need, start and end time of the event.  We wrote python code that looked for the jump in signal value and return to baseline.  The trouble is that each machine has a different baseline value.

Are there simple tools/techniques that can determine from recent context the end of a pronounced signal?  My preference is python but I'm happy to use anything that works.  Thanks in advance for any feedback.

This sounds like a job for a simple high-pass filter, followed by a simple threshold, to me.

1/3 billion records isn't much data if you're used to radio signal sampling rates, so computationally, this should be most definitely doable.

So, with scipy.signal there's a reasonable signal processing library for Python, but I'm from the GNU Radio world, where you can use Python to plug together a signal processing flow graph.

In your case, it'd look like this:

Vector Source -+> high pass filter --> threshold --> [custom block]

which  you can click together graphically in GNU Radio companion and let it generate Python, or you can directly write in Python.

The [custom block] is from the top of my head something you'd have to write yourself. You can do that in Python or C++, but Python will be slower to run. I'd still do it in Python for a first prototype:

When the first input is 0, it simply consumes the the input (effectively dropping that data). When it's 1, save the number of the first 1-sample in a row. Then wait until the 1s stop and safe locally something

self._list_of_signals.append( (first_item_number, last_item_number ) )

I think you get the idea.

You feed in the data into the vector source in your python framework, e.g. with the data from your database, and then let the flow graph run, i.e. process all the data it has.

When it's done running, you ask your custom block to give you the list of timestamps. With their first sample number times the sampling rate, you get a time stamp, and then you just put (start, stop) into your database of "trimmed" data.


Generating audio clip by superposing two files I'm currently trying to generate a set of audio files where some audio events (e.g. dog bark) are immersed in a background audio scene (e.g. crowd noise).

I understand that if the event would have happened in the original scene, then there would also be some components due to the interaction between the sound and the environment. In other words, this should be, from my understanding, a convolutive mixture.

My question is, if I approximate the mixture by summing the two signals, thus doing an instantaneous mixture, am I losing a crucial part of the audio signal? In which circumstances is this a reasonable assumption, if any?
My final objective is to train a classifier with the resulting audio.

We typically model audio to be a pretty linear thing, so yes, just add – you can multiply each input with a constant to weigh them differently.


How can I know the type of filter from its cutoff frequency? So I have the following equation
$$
y[n]=0.9y[n-1]+0.1x[n]
$$

We can find easily find the transfer function, which is

$$
H(z)=\frac{0.1}{1-0.9z^{-1}}
$$

and from that, the frequency response, which is

$$
H(e^{j \omega})=\frac{0.1}{1-0.9 e^{-j \omega}}
$$

To find the cutoff frequency we just make the magnitude of the above equal to $\frac{1}{\sqrt{2}}$. Then, $\omega_0 = 0.0675\Rightarrow f_0 = 0.0107$.

Okay, so how can we say that given  $f_0 = f_c$ the filter is a LPF or HPF?

Types of filters, such as LPF, HPF, BPF or BSF, are not described by their cutoff frequencies $f_c$ but through their passband or stopband frequency $\omega$ region.

For example, in your case of the filter $H(z) = 1 / (1 - 0.9 z^{-1})$, the passband is in the low frequencies (where the magnitude of the frequency response $H(e^{j\omega})$ is greater than $1/\sqrt{2}$. And by your definition that region corresponds to $\omega \in (-f_c , f_c)$ which indicates a lowpass range of frequencies.

Other filters like HPF and BPF follow in a similar manner. Note that specific filters such as comb or notch types are defined in a different style than this.


Performance of adaptive filters Can somebody please provide an intuitive answer or reference for the following questions?

Q1: Dependence of estimation performance on number of data points -- I could not find any information whether the estimation performance of Adaptive filters such as Least Mean Square (LMS), Constant Modulus Algorithm (CMA) and Kalman Filters depend on the number of data points or not. Is there any information whether the mean square error between the actual and estimated parameters reduces with the increase in the number of data points or not?

Q2: Dependence of convergence on number of data points -- For instance, information such as if convergence of these adaptive filters (or in general) depends on the number of data points i.e., if these require a large number of data points to have good estimation performance.

Q1: Dependence of estimation performance on number of data points

Since LMS and RLS are adaptive filters, their estimation performance improves as the number of their iterations increase. Hence more data points will make their outputs closer to the expected performance, until the convergence is achieved (this requires either a WSS data or a slowly changing, statistically,  nonstationary data, so that filter can track it's statistical character). Once the filter is operating in convergence conditions, then there won't be any improvements in its output, providing what is called as the (minimum) steady-state error. This steady-state error depends on a number of factors but not on the data length.

But until the convergence is achieved, the estimation error decreases as more data points processed. However most typical adaptive filter users would be interested in its steady state (convergent) results rather than the transient response.

Note that increased data points can either come from a longer observation or from higher sampling rate, the results of which can be different. For example for a Kalman filter (in its extended mode) with mechanical applications, an increase in sampling rate can actually improve the steady state error as well.

  Q2: Dependence of convergence on number of data points

For both LMS, RLS and Kalman filters, convergence primarily depends on the number of data points being processed; i.e., number of iterations. However this convergence rate is different for different filter types, reflecting their complexity and/or sophistication. The simple LMS filter has the slowest rate of convergence (roughly 10 to 20 times the filter tap weight length) whereas the RLS and Kalman filters display a convergence rate of roughly $2$ times the filter length (Haykin_Adaptive Filter Theory) for WSS inputs.


Why Wiener filter not IIR I would like to limit discussion to discrete time causal Wiener filter.

By definition it is an FIR filter, and it also has optimality under Mean Square Error criterion within the class of linear filters.

How do we know that there isn't a better IIR filter (something with a feedback loop), which is also linear?

You're right, there usually is a better IIR filter (if you have enough data). The discrete-time Wiener filter is not "by definition" FIR. It is common to constrain the filter to the FIR case because it's often more straightforward to implement, and because such a filter can be made adaptive more easily. Also, in practice you often want to consider only a finite data window, which makes the optimum Wiener filter FIR.

However, you don't need to impose the constraint that the filter length be finite. The most general causal discrete-time Wiener filter for data available over the infinite past is an IIR filter.

A derivation of the FIR and the IIR cases can be found (among others) in Digital Signal Processing by Proakis and Manolakis.

EDIT: If we only have a finite amount of data available, the optimum Wiener filter has a finite impulse response. Assume we have $N$ data points (one current and $N-1$ consecutive past data points), then the best we can do with a linear time-invariant filter is to linearly combine those $N$ data points, by multiplying each data point by a coefficient and adding them up:

$$y[n]=\sum_{k=0}^{N-1}h[k]x[n-k]\tag{1}$$

The coefficients $h[k]$ are optimized such that the mean square error between the filter output $y[n]$ and some desired signal is minimized. The $N$ coefficients $h[k]$ are the impulse response of the corresponding FIR Wiener filter.


Is there a generalized method to get the input with the given output and the impulse response? $y(t)=y(t+12), y(t) = x(t) \ast h(t)$

The continuous time signal output $y(t)$ is a periodic square wave, 50% duty cycle pulse.

The impulse response is a box function.($A = 1, T = 2$)

By using Fourier series and Fourier transform, I thought I can guess the input $x(t)$, assuming the input is also periodic with the period 12 as well.

$$
w_0=\frac{2\pi}{12}\\
x(t) = x(t+12) = \sum_{k=-\infty}^{\infty} a_k e^{jkw_0t},
a_0 = \frac{1}{T} \int_{0}^{T} x(t)dt,
a_k = \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} x(t) e^{-jkw_0t}dt \\
y(t) = y(t+12) = \sum_{k=-\infty}^{\infty} b_k e^{jkw_0t},
b_0 = \frac{1}{T} \int_{0}^{T} y(t)dt,
b_k = \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} y(t) e^{-jkw_0t}dt \\
$$

By using the Fourier transform,
$$
X(jw) = 2\pi \sum_{k=-\infty}^{\infty} a_k \delta (w-kw_0)\\
Y(jw) = 2\pi \sum_{k=-\infty}^{\infty} b_k \delta (w-kw_0)\\
H(jw)=2sinc(w) \\
X(jw) = Y(jw)/H(jw)
$$

I was stuck at this point; I thought I can write an equation about these coefficients, but the sinc function was a hassle for me.

How should I write an equation which is related to $a_k, b_k$?
Also, in order to guess the input, do I have to use the Fourier Transform like this?
Or is there a characteristic of the convolution for periodic/box functions?

I wasn't sure whether my approach was correct to find inputs with the given output and the impulse response. Thanks for your time.

The output signal of an LTI system is periodic only if its input signal is periodic as well, your reasoning was right. Now, a signal periodic with a certain period $T_p$ has a line spectrum, i.e., it only contains frequency contents at $n f_0$, $n \in \mathbb{Z}$ where $f_0 = \frac{1}{T_p}$. What it basically means is that your input only excites these frequencies so that to obtain the input output relation, it only matters how the LTI system responds to these particular frequencies. You can therefore consider the impulse response of your LTI system periodic as well, it changes nothing.

That said, you can express the input-output relation entirely using Fourier series. If your input signal $x(t)$ has Fourier coefficients $C_{x}[\mu]$, your output signal $y(t)$ has Fourier coefficients $C_{y}[\mu]$ and your system is characterized by an impulse response $g(t)$ with Fourier coefficients $C_g[\mu]$ then, the input output relation is

$$y(t) = x(t) * g(t) \quad \Leftrightarrow \quad C_y[\mu] = C_g[\mu] \cdot C_x[\mu],$$

where * represents circular convolution. Therefore, if you are given $y(t)$ and $g(t)$, you can find $x(t)$ via:

Obtaining the Fourier coefficients $C_y[\mu]$ of $y(t)$
Obtaining the Fourier coefficients $C_g[\mu]$ of the (periodified) $g(t)$
Obtaining the Fourier coefficients $C_x[\mu] = C_y[\mu]/C_g[\mu]$
Finding $x(t)$ via its Fourier series expansion using $C_x[\mu]$.

Let me elaborate shortly on the difference between Fourier transforms and Fourier series: The Fourier transform takes a continuous time signal $g(t)$ and returns its amplitude density $G(\omega)$ for every $\omega$. You find it via $$G(\omega) = \int_{-\infty}^\infty g(t) {\rm e}^{-\jmath \omega t} {\rm d}t.$$ On the other hand, the Fourier series is defined for periodic functions only. It takes a periodic time signal and returns coefficients belonging to the integer multiples of its base frequency. You find the Fourier coefficients via $$C_g[\mu] = \frac{1}{t_p} \int_{t_p} g(t) {\rm e}^{-\jmath 2\pi \mu t/t_p} {\rm d}t.$$ if $g(t)$ is $t_p$ periodic. These two definitions may look very similar. Besides normalization it may seem that $C_g[\mu]$ is the same as $G(\omega)$ at $\omega = \mu \cdot \frac{2\pi}{t_p}$. However, this is only true if $g(t)$ is either actually periodic or at least time-limited to the interval $t_p$ you are considering (which is true in your case!) since in this case, periodifying it does not change it at all. If it is not limited to $t_p$, the Fourier series integral limited to the interval $t_p$ would truncate the function which would lead to wrong results. Therefore, one would have to periodify $g(t)$ first in this case to get correct results.


relationship between energy and sampling rate Excuse my silly question, but i really want to know does changing the sampling rate affects the energy (bandwidth) of a signal? therefore improves cross correlation output?

You must observe the nyquist frequency when sampling a signal.
In order to sample a signal without introducing artifacts, you must first filter out everything in the signal that has a frequency higher than half of your sampling rate.  If you sample a 1000Hz, then you must first filter out everything above 500Hz.

Yes, changing the sampling rate can change the bandwidth of digital representation of the signal.  A lower sampling rate means a smaller bandwidth.

This will only improve crosscorrelation if the noise is all high frequency stuff, and the signal you are trying to detect is in the lower frequencies.  You could get the same effect, however, by using a low pass filter  on the sampled data.  That would also remove the high frequencies and make the correlation clearer.


BIBO Stability for system with no poles I have two questions regarding systems with no poles:

Why does a system with no poles have a finite support?
Why if the system has a finite support it means that it is BIBO
stable?

Why does a system with no poles have a finite support?

If a system doesn't have finite poles, then its transfer function is of the form:

$$H(z) = \frac{Y(z)}{X(z)} = a_Nz^N+a_{N-1}z^{N-1}+...+a_1z+a_0$$

So if you go back to the time domain:

$$y[n]=a_Nx[n+N]+a_{N-1}x[n+N-1]+...+a_1x[n+1]+a_0x[n]$$

Note that if $x[n]=\delta[n]$, then $y[n]=h[n]$, the impulse response of the filter.

$$h[n]=a_N\delta[n+N]+a_{N-1}\delta[n+N-1]+...+a_1\delta[n+1]+a_0\delta[n]$$

It's easy to see that the equation above will return a $h[n]$ that is non-zero for a limited range of values of $n$.

This means that the filter is FIR. Therefore, it has a finite support and it is BIBO stable, leading to your second question.

  Why if the system has a finite support it means that it is BIBO stable?

For BIBO stability of an LTI discrete-time system, one must check that the following condition is satisfied:

$$\sum_{n=-\infty}^{\infty}|h[n]|<\infty$$

If $h[n]$ is the impulse response of a FIR filter, then it has some length $L$:

$$\sum_{n=-\infty}^{\infty}|h[n]|=\sum_{n=K}^{L+K-1}|h[n]|$$

Due to this summation being finite, it converges. Ergo, any FIR filter is BIBO stable.


Compute step response from impulse response of continuous-time LTI system The impulse response of the LTI system is

$$h(t)=e^{-4t} u(t)$$

The expression for the step response is

$$\frac14 \left(1-e^{-4t}\right)u(t)$$

My question is how $u(t)$ appears in the answer.

Since your question is not about how to obtain the answer, but specifically about the reason why the unit step $u(t)$ appears in it, I assume that you know that the response to a step input is computed by convolving the impulse response with a step:

$$y(t)=\int_{-\infty}^{\infty}h(\tau)u(t-\tau)d\tau=\int_{-\infty}^{t}h(\tau)d\tau\tag{1}$$

Since $h(t)$ is causal, i.e., $h(t)=0$ for $t<0$, the lower integration limit can be changed to $0$. This implies that if in $(1)$ the upper integration limit $t$ is less than $0$, the result must be zero. Consequently, the response to a step input can be written as

$$y(t)=u(t)\cdot\int_0^{t}h(\tau)d\tau\tag{2}$$

Evaluating $(2)$ gives you the final answer, which includes the step function $u(t)$.

EDIT: Triggered by Dilip Sarwate's comment I add some extra explanation. Note that the step function in $(2)$ is actually redundant because for $t<0$ the integral in $(2)$ is zero because for $t<0$ we have

$$\int_0^{t}h(\tau)d\tau=-\int_{-|t|}^0h(\tau)d\tau=0,\quad t<0$$

because $h(t)=0$ for $t<0$. However, if $h(t)=f(t)u(t)$ with some function $f(t)$ that is not zero for $t<0$ (in the given example we have $f(t)=e^{-4t}$), we cannot just replace the lower integration limit in $(1)$ by $0$ and leave out the step function:

$$\int_{-\infty}^th(\tau)d\tau=\int_{-\infty}^tf(\tau)u(\tau)d\tau=\int_{0}^tf(\tau)u(\tau)d\tau\neq \int_{0}^tf(\tau)d\tau$$

because now we generally have for $t<0$

$$\int_{0}^tf(\tau)d\tau=-\int_{-|t|}^0f(\tau)d\tau\neq 0$$

The correct way to evaluate the integral for $h(t)=f(t)u(t)$ is

$$\int_{-\infty}^th(\tau)d\tau=u(t)\cdot\int_0^tf(\tau)d\tau$$


What does the operator 'Range Detector' do in this block diagram? How can it be modeled? I'm trying to implement a 'Compressor/Expander/Noise gate' for automatic gain control of audio signals. The block diagram of the system is shown below:

I don't understand how the 'Range Detector' selects the appropriate output, given that for some values of xRMS(n), both the output of the expander and noise gate are non-zero.

According to the source that your block diagram is taken from [Udo Zölzer - Digitale Audiosignalverarbeitung], the Range Detector is an arbiter whose decision is not based on the output of the three nodes, but rather on the region of operation of the static characteristic that your calculated $x_{RMS}(n)$ is found in:

So, if $x_{RMS}(n)$ is below $NT$, the range detector chooses the lower node, if it is above $CT$, it chooses the upper node...

As you already noticed, that block diagram does have a few shortcomings when it comes to implementation:

It does not take into account edge cases like $|ET| > |CT|$
There are more regions of operation than branches to process the RMS value. Even in the case where your RMS value is above expander threshold but below compression threshold you would be required to perform a set of logarithmic, multiplication, comparison, addition and exponential operations.
But it would be computationally more efficient to calculate whether $2 ^ {2 * ET} < x_{RMS}(n) < 2 ^{2 * CT}$ evaluates to true and if yes, you can set $f(n) = 1$ immediately without having to calculate $log_2(x_{RMS}(n))$ and so forth.
The lowest branch is not actually a noisegate but another expander with a (presumably) steeper slope.  Zoelzer defines the slope as $S = 1 - \frac{1}{R}$ while a Noisegate's ratio is defined as $R = 0$.
As you can see, the slope would thus evaluate to $NS = -\inf$ and the calculation of the whole lower branch becomes rather pointless.
Again, it would be more efficient to set $f(n)$ directly to $0$ in the case $x_{RMS} < 2 ^{2 * NT}$ and save yourself a few pointless calculations.

As a matter of fact, Zoelzer seems to agree with the last point and provided a sligtly different block diagram in the book DAFX:

Please excuse the poor image quality. I highly recommend you get this book as well. It is a great addition to Digital Audio Signal Processing and provides Matlab demo code for most topics, dynamics processing included.


Edge map based on a Haar Wavelet Transform I have been implementing the paper  Blur Detection for Digital Images Using Wavelet Transform and was asking myself how the following formula could reconstruct the edges given a Haar Wavelet transformed image :

$$ \sqrt{LH_i^2 + HL_i^2 + HH_i^2} $$

This formula is on page 2, Algorithm 1, step 2.

I looked at the reference given, but couldn't find where this was explained.

One single level of a standard separable 2-channel wavelet transform, denoted by $i$, uses a low-pass $l$ and a high-pass $h$  filters (followed by downsampling). Traditionally, one applies $l$ and $g$ on the rows of the image, putting the downsampled low-passed coefficients on a left-half, and the downsampled high-passed coefficients on a left-half. Then the resulting image is processed in the same way column-wise.

This results in a set of four separable 2D filters:

$l$ on rows, $l$ on columns (LL)
$l$ on rows, $h$ on columns (LH)
$h$ on rows, $l$ on columns (HL)
$h$ on rows, $h$ on columns (HH)

producing the arrangement below. This is iterated (levels) for on the image, then on the low-pass/low-pass images (LL)

Among those filters, three (LH, HL, HH) can be considered as edge detectors ($h$ is similar to a 2-point gradient) in different directions: more horizontal (HL), more vertical (LH), somehow diagonal (HH).

The corresponding coefficients, as least their combined energy, can be interpreted as a measure of edge strength, a map computed for each level.

Remember that for a continuous 2D field modeling an image $f(x,y)$, the image gradient

$$\nabla f =\left[\frac{\partial f(x,y)}{\partial x} ,\frac{\partial f(x,y)}{\partial y}  \right]^T$$

represents the local directional change in image intensities, whose norm is a measure of gradient "magnitude". You can compute several norms, the most common being the Euclidean, the Taxicab or the max norms.

The Haar filter $[1\,-1]$ is one possible discretization for the 1D gradient, the three 2D Haar wavelets compute an horizontal, a vertical and a somewhat diagonal gradient, hence their combined energy (Euclidean norm) is an estimate of the true gradient magnitude.


Is there an efficient algorithm to extract low frequency signal from audio with high frequency carrier? I need to process an audio signal in an Android app. As shown in the two examples below, there is a 16kHz carrier signal and a lower frequency signal of varying frequency (0-4kHz). Does anyone know of an efficient algorithm to extract the lower frequecy signal in an Androd app?

http://mc.weatherflow.com.s3.amazonaws.com/media/IMG_03102014_135912.png

If you know your carrier is always exactly 16kHz (within a few percent or so), and you have sampled the audio at 48kHz, you can try a simple 3-tap FIR filter approach:

$$Smoothed(k) = \frac{1}{3}( Sig(k) + Sig(k-1) + Sig(k-2) )$$

This assumes your original signal is called $Sig()$, and it computes a new signal, $Smoothed()$.

The method exploits the fact that three sample periods equates to one entire cycle of a  16kHz sine wave, so an 16kHz component in the input signal will be cancelled (since summing all the samples of a full sine wave cycle will result in zero).

The resulting $Smoothed()$ signal should look like you original plots, with the furry 16kHz components gone.

If your sample rate is not 48kHz (let's call it $F_s$), or if your interfering sine wave is at some other frequency (let's call it $F_i$), you can craft the 3-tap filter as follows:

$$Smoothed(k) = \frac{1}{2-2\cos(\omega)}( Sig(k) - 2\cos(\omega) Sig(k-1) + Sig(k-2) )$$

where $\omega = 2\pi\frac{F_i}{F_s}$ (note that, if $F_i=16000$ and $F_s=48000$, we get $\omega=\frac{2\pi}{3}$ and hence $\cos(\omega)=-\frac{1}{2}$. If you substitute this into the second equation, you will simply get the first equation again.

This method is simply creating a $2^{nd}$ order FIR filter with zeros at $\pm F_i$, and with the gain at DC normalised to unity.


"Normalizing" my D4 wavelet transformation at each step reduces final image quality Original image:

(Images included are .png images, so no additional distortion was added in save/upload for viewing)
I have used the D4 transform from page 20 of "Ripples in mathematics", which is basically these 5 steps:
Forward d4:
c1 = √3 / 4.0 ;
c2 = (√3 - 2) / 4.0 ;
s[ IEVEN ] += √3 * s[ IODD ] ;
s[ IODD ] -= c1*s[ IEVEN ] + c2*s[ IPREVEVEN ] ;
s[ IEVEN ] -= s[ INEXTODD ] ;
s[ IEVEN ] *= ( √3 - 1 ) / √2 ;
s[ IODD ] *= ( √3 + 1 ) / √2 ;

The inverse:
c1 = √3 / 4.0 ;
c2 = (√3 - 2) / 4.0 ;
s[ IODD ] *= ( √3 - 1 ) / √2 ;
s[ IEVEN ] *= ( √3 + 1 ) / √2  ;
s[ IEVEN ] += s[INEXTODD] ;
s[ IODD ] += c1*s[ IEVEN ] + c2*s[IPREVEVEN] ;
s[ IEVEN ] -= √3 * s[ IODD ] ;

I'm compiling and running this using double precision values from C++.  I run this on the rows of the image, then the columns.  I use a crude filtration algorithm to remove the lowest 90% of the difference coefficients in the image.
The filtration algorithm is:

Run through entire transformed image (as a set of numbers)
Find the largest difference coefficient (maxVal) (in entire 2d image)
Choose minValToSurvive as 1% of maxVal.
If a difference coefficient has a magnitude less than minValToSurvive, it is zeroed.

Here's my problem is.  When I remove only 83% of the lowest difference coefficients from the image (minValToSurvive=0.01*maxVal), you get this:
normalized

If I remove the normalization steps:
s[ IEVEN ] *= ( √3 - 1 ) / √2 ; // REMOVE
s[ IODD ] *= ( √3 + 1 ) / √2 ;

(in both the fwd and reverse transforms), the result after removing 90% of the components is much better (much less noise)

So I can think of 1 of 2 problems:

Normalizing the image by the ( √3 - 1 ) / √2 factors is killing precision
I'm not filtering correctly

Or am I wrong?  If I'm filtering (removing insignificant components) incorrectly, what is a better way to filter?  If it's the floating point precision, then should I not normalize the transform at every step?

The right answer is you have to perform each of the Update/Predict steps on the input signal completely before doing the next Update/Predict.  What I was doing was walking through the signal, and performing each Update/Predict as I go.

On page 158 of "Ripples", there is a reference implementation.

// s is the signal
#define IEVEN (2*j)
#define IODD (2*j + i)
for( int i = 1 ; i <= n/2 ; i *= 2 )
{
  for( int j = 0 ; j <= n/2 - i ; j += i ) // Must do this Predict step COMPLETLEY
    s[ IEVEN ] += √3 * s[ IODD ] ;

  for( int j = 0 ; j <= n/2 - i ; j += i ) // Then this one..
  {
    int prevEvenIndex = IPREVEVEN ;
    s[ IODD ] -= d4c1*s[ IEVEN ] + d4c2*SAFE_PREV(s,prevEvenIndex) ;
  }

  for( int j = 0 ; j <= n/2 - i ; j += i )
  {
    int nextOddIndex = INEXTODD ;
    s[ IEVEN ] -= SAFE_NEXT(s,nextOddIndex) ;
  }

  for( int j = 0 ; j <= n/2 - i ; j += i )
  {
    s[ IEVEN ] *= d4normEvens ;
    s[ IODD ] *= d4normOdds ;
  }
}

The 98% 0's D4 transform:


use chirp z transform for spectral resolution let us suppose that we some signal $x(t)$,which consist of sinusoidal   components and white noise,i would like to know how to use    chirp z transform for improving spectral resolution?i found   this article about  chirp z transform

http://prod.sandia.gov/techlib/access-control.cgi/2005/057084.pdf

also  similar  article on official matlab side

http://www.mathworks.com/help/signal/ref/czt.html

but  as i know always in spectral estimation methods,some  knowledge we should know how process is going,or  some other things,let us suppose that we know sampling frequency lat say $Fs=100$,and   also we know that it consist as i mentioned by periodic components(sin or cosine),but i dont  know number of components,neither frequency or neither phases(+neither amplitudes),how can i apply chirp z transform to this signal?we have just some signal samples  with size $N=294$ ,please help me how to  choose correct components so that  get  effectively spectral structure and if let say one frequency is $f_1=13$ and second $f_2=13.5$,this transform be able to separate this  frequencies from each other,thanks

If you have no prior knowledge about the approximate locations of the frequencies, the Chirp Z-transform is of no immediate use to you. The Chirp Z-transform functions like a magnifying glass, so you need to know where you want to look and the Chirp Z-transform will show you the details. I would suggest you use an FFT to get an idea where the frequencies are, and if you need a very high resolution in a certain area of the spectrum, then the Chirp Z-transform can be useful. However, if computational complexity and memory are no issue, you can simply use a long FFT to achieve the same. If you want to know more about the Chirp Z-transform I would recommend this paper to you.


Why discrete wavelet transforms use sampling rate 2 and need signal length to be a power of 2? I know Fourier transforms but new to wavelet transforms. I can understand Haar transform needs signal length a power of 2, since the filters have 2 taps and down-sampling and up-sampling in the pyramid algorithm form binary tree structures, so it is efficient to implement if signal length is a power of 2.
But higher-order transforms (e.g. DB2) or other kinds also use sampling rate of 2 and need signal length a power of 2. Is it also for efficient implementation? Can the sampling rate be other small prime numbers 3, 5 like in FFT?

I can't comment because of reputation but if i'm not mistaken, the need for DWT to be supplied an input that is length power of 2 is independent of it's subsequent levels being downsampled by a power of two. The input does not have to be a power of two but you'll have to deal with problems at the edges of the signal then.
Here is a good link for reference that addresses padding of a signal for wavelets:
https://paos.colorado.edu/research/wavelets/faq.html#padding
Basically states that to prevent problems at the edges of the tranform, the inputs are usually padded.
I've also seen places mention other decimation schemes but it begs the question as to why they would be necessary?
For every level of decomposition in a DWT, the signal is downsampled by 2 because the subsequent level of decomposition is limited to half the bandwidth of the previous level. Allowing for half of the samples to be discarded.


What are poles and zeros? The concept of poles and zeros in filters was introduced to me during our lab (our lab isn't sync with our lecture) through the pole-zero plot generation of filters in MATLAB. I didn't get its importance until it was mentioned in our lecture that poles make the equation infinity (or its denominator zero) while the zeros make the roots in the numerator zero. I still can't quite grasp what poles and zeros are. Can someone give an an intuitive and simple explanation on what poles and zeroes are? Much appreciated.

Take the equation b/(x-c) with b non-zero.  The ratio goes to infinity as x approaches c.  So c is the location of a pole (something tall and pointy in a graph).

Take the equation (x-b)/c with c non-zero.  The ratio goes to zero as x gets closer to b.  So b is the location of something commonly called "a zero".

You can not only do this with scalar x, but with complex x, thus the domain of the poles and zeros will be on a (complex) plane instead of on a line.

If the ratio represents something about the response of a filter, it might say the filters output is at or near zero when the input is at or near a response "zero".  And bad things might start to happen when x gets close to a pole (power supply starts smoking when asked to supply infinite amps, math operations produce NaNs or fixed-point overflows, etc.)


Dimension of Matrix in MIMO with multi-path In a MIMO system, the received signal is supposed to be $y = Hx + n$, where $H$ is the channel of dimension $N_t \times N_r$, $x$ is the transmitted signal with dimension of $N_t\times 1$ , and $n$ is the noise, finally the resulted signal $y$ represents the received signal of dimension of $N_r\times1$.

Assuming that I need to use channel $H$ as multi-path channel whose number of taps is $e.g.$ 5 .. The question is:

Will the dimension of channel $H$ will be $5N_t\times N_r$?
If so, how can we perform the matrix multiplication of $Hx$ in  $y = Hx + n$, so the received signal $y$ will be of dimension $5N_r\times1$?

The short answer is: if you transmit one block of $N$ symbols (no data before it for at least $L$ symbol times) over a frequency-selective channel of length $L$, then the channel matrix will be of dimension $NN_TN_R\times N$.

The long answer: start with $N_T=N_R=1$. Then the $n^{\text{th}}$ received sample can be written as

$$y_n=\sum_{l=0}^Lh_lx_{n-l}+z_n$$

If you arrange these samples in a vector-matrix form you get

$$\mathbf{y}_{N\times 1}=\mathbf{H}_{N\times N}\mathbf{x}_{N\times 1} + \mathbf{z}_{N\times 1}$$

Now let $N_T=1$ while the number of receive antennas $N_R$, then the $n^{\text{th}}$ received sample over $m^{\text{th}}$ receive antenna is

$$y_n^{(m)}=\sum_{l=0}^Lh_l^{(m)}x_{n-l}+z^{(m)}_n$$

for $m=1,\,2,\,\ldots,\,N_R$. If you arrange these sample in a vector-Matrix form, you get

$$\mathbf{y}_{NN_R\times 1}=\mathbf{H}_{NN_R\times N}\mathbf{x}_{N\times 1} + \mathbf{z}_{NN_R\times 1}$$

Similarly, you can do the general case for $N_T$ transmit and $N_R$ receive antennas.

As you can see, to make this work on block-by-block basis, you need to padding $L$ zeros before each block to eliminate the inter-block interference. Alternatively, you can use cyclic-prefixed (CP)-OFDM, which is more efficient at the receiver side.


Identifying origin in convolution table I am taking the convolution of $x(n)=\left\{2,1,-1,-2,3 \right\}$ with $n=0$ at the third position with $h(n)=\left\{1,2,0,3 \right\}$ with $n=0$ at the second position. The answer is $y(n)=\left\{2,5,1,-10, -10, -3, 6, -9 \right\}$ with the $n=0$ at the fourth position. I studied convolution more than a decade ago and I used a table to do the convolution sum as follows:

There used to be a simple rule to obtain the location of the $n=0$ position of the output $y(n)$ from the table itself or perhaps from $x(n)$ and $y(n)$. But alas I cannot remember it. Grateful for any suggestions.

Given a signal $x[n]$ with the nonzero range of $N_1 \leq n \leq N_2$ and the signal $h[n]$ with the nonzero range $M_1 \leq n \leq M_2$, then the nonzero range of their comvolution $y[n] = x[n] \star h[n]$ is $N_1+M_1 \leq n \leq N_2 + M_2$.

Applying this to your case yields the following:
For the signal $x[n]$ with $n=0$ at third position the nonzero range is $-2 \leq n \leq 2$ and for $h[n]$ with $n=0$ at the second position the nonzero range is $-1 \leq n \leq 2$ which therefore yields the nonzero range of their convolution $y[n]$ is:
$$ -3 \leq n \leq 4 $$

which means that the $n=0$ for $y[n]$ is at the fourth position.


Image Processing on Raspberry Pi What are my options for Open Sourced Image Processing functionality on a Raspberry Pi?

So, I am referring to the old-school processing, using operators like open/close/erosion/dilation/threshold/convolution etc.

I am not looking for neural network based processing.

Bonus points if it has a Python interface, but C will work as well.

In the 1990s, in University I was using Scil Image for this stuff, and I wonder what the modern equivalent for it is?

Python has a lot of import libraries for image processing, and they are running on platforms that support Linux, such as Raspberry Pi.

Pi runs with ARM based processors, and every limitation of that CPU (especially on SIMD floating point ops) is also a limitation for image processing.

Pi also comes with a dedicated hardware GPU which can be used to accelerate certain image processing (computer graphics) functionality (in addition to its image/video encode and decode functionality)

Pi also comes with built in camera interface that would enable raw processing of camera data, using available Python libraries.


Design a comb notch filter when sampling rate is not an integer multiplicative of the cutoff frequency Let's say I have a signal sampled at 976 Hz and I want to remove network interference noise, i.e. 50 Hz and all the harmonics. Is it possible theoretically to design a comb filter to do that?

I have looked at how to do this in MATLAB but the problem is that the order of the comb filter -which defines the fundamental frequency to be filtered out- must equal fs/fo, and of course this needs to be an integer.

Any help would be greatly appreciated.

You can design that filter manually without problems. Matlab just uses a very simplistic approach to comb filtering with a delay line.

In order to keeps things as simple as possible I would recommend you use a series of notch filters to remove each partial of your harmonic noise separately. That also gives you more control over how much of each harmonic really has to be removed.

A simple notch filter is one with unit gain at DC, i.e. $H(1) = 1$, and the zeros on the unit circle at the desired notch frequency, positive and negative. The poles have to be inside the unit circle, close to the zeros to cancel their effect farther away from the notch frequency. We can design this directly in the z-domain.

$$H(z) = \frac{(z-\exp(i\omega))(z-\exp(-i\omega)}{(z-r\exp(i\omega))(z-r\exp(-i\omega))}\times\frac{(1-r\exp(i\omega))(1-r\exp(-i\omega)}{(1-\exp(i\omega))(1-\exp(-i\omega))}$$

Here $\omega$ is the normalized frequency, meaning $\omega=2 \pi \frac{f_0}{f_s}$ where $f_0$ is your notch frequency and $f_s$ the sampling frequency. The parameter $r\in]0,1[$ controls the width of the notch. The second fraction makes sure we get a unit response at $\omega=0$.

To get the coefficients for a recursive filter implementation, we have to simplify the transfer function and cancel $z^2$ to turn it into a rational function of $z^{-1}$. The result is
$$H(z) = \frac{\left(1+r^2-2r\cos(\omega)\right)\left(z^{-2}+1-2z^{-1}\cos(\omega)\right)\csc(\omega/2)^2}{4\left(r^2 z^{-2}+1-2rz^{-1}\cos(\omega) \right)}$$

So you can read off your filter coefficients to be

$$ A[n]=\left(1,-2r\cos\omega,r^2 \right)$$
$$ B[n]=\left(1,-2\cos\omega,1 \right)\cdot\frac{1}{4}\left(1+r^2-2r\cos\omega\right)\cdot\csc(\frac{\omega}{2})^2$$

The last open question is how exactly the parameter $r$ controls the notch behavior. Near the zero the transfer function on the unit circle locally behaves like

$$f(x) = \frac{x}{x-i(1-r)}$$

around $x=0$. For $r=1$ we get a singularity at $x=0$ and an ill-defined transfer function. For values $r<1$ but close to $1$ we can see that for both $x\to-\infty$ and $x\to\infty$ the transfer function becomes $1$. The interesting question is, how far does the effective influence of the pole reach on the x-axis?

The square magnitude of $f(x)$ is $$|f(x)|^2=f(x)\cdot f(x)^*=\frac{x^2}{x^2+(1-r)^2}$$
which is always real and positive. We can easily equate it to a threshold and solve for $x$.The most common threshold for filter cutoffs in signal processing is $1/\sqrt{2}$, and if we solve $|f(\Delta x/2)|^2==1/\sqrt{2}$ we get the effective bandwidth of the notch to be
$$\Delta x = 2(1-r)\sqrt{1+\sqrt{2}}$$

and we can use the equivalence of our problem to the original pole/zero placement to conclude

$$ r = 1-\frac{\Delta \omega}{2\sqrt{1+\sqrt{2}}} $$

with the natural frequency bandwidth $\Delta \omega$. If you want this as a true frequency then again $\Delta \omega=2\pi\frac{\Delta f}{f_s}$

Note that this bandwidth calculation used a local approximation, so it is strictly only true for $(1-r)<<1$, or small bandwidths compared to the sampling frequency.

So this gives you a notch filter for one single harmonic of your harmonic noise. Just apply one filter at a time, one after the other, with the parameters matched to the harmonic you would like to cancel, up to fs/2. The matlab filter() function readily takes the coefficients A[n] and B[n] given above.

An attempt to implement this in scilab seems to work nicely.  Plot of before and after spectrum:

Code:

// 12740
fs = 1000;
f = 101.89798;
omega = 2*%pi*f/fs;
delta_omega = omega/2;
N= 1000;
t = [0:N];
f = [0:N]/N*fs-fs/2

r = 1 - delta_omega/(2*sqrt(1+sqrt(2)));

num = [1 -2*cos(omega) 1]*1/4*(1+r.^2-2*r*cos(omega))*csc(omega/2).^2;
den = [ 1 -2*r*cos(omega) r.^2];

X = sin(omega/2*t) + sin(omega*t) + sin(2*omega*t);
Y = filter(num,den,X);

figure(1)
clf;
plot(f,log(fftshift(abs(fft(X)))+0.01))
plot(f,log(fftshift(abs(fft(Y)))+0.01),'k:')


How can I get the power of a specific frequency band after FFT? I need help with my matlab code. I have recorded EEG data and now I want to analyse it in matlab. I want to find the mean power of specific frequency ranges (theta: 4-8 Hz,alpha: 8-12Hz etc.). For that I performed a fft. Next I used the interp1 function to find the amplitudes at specific frequency ranges. But I am not sure if I have done it correctly...It would be great if someone can look at the code.

Thank you for your help!

%p2 is the data in the time domain
% Variables

fftlength = 1024;                                % Set FFT length
stepsize = 32;                                   % Set step size
samples = fftlength:stepsize:length(p2);    % create an index array
c = 1:fftlength;

%Specify hanning window

hanning = [1:fftlength]';
hanning_in = 2* pi() * (hanning - (fftlength+1)/2)/(fftlength+1);

%rescaled x-axis to match sample length?

hanning1 = (sin(hanning_in)./hanning_in).^2;

f=[128/fftlength:128/fftlength:128];   % frequency index for the spectral
array

bp_p2=cell(1,1);

channel=1;

for kk = 1:length(samples)

k = samples(kk);

spectrum = fft(p2(k-fftlength+1:k,:) .* hanning1)/fftlength; % apply window
to filtered data

spectrum2 =(sqrt(spectrum .* conj(spectrum)));  % get magnitude uV

% theta (4-8),alpha (8-12),beta low(12-16), beta high(16-25)

theta=mean(interp1(f,spectrum2(:,channel),[4:f:8]));
alpha=mean(interp1(f,spectrum2(:,channel),[8:f:12]));
betalow=mean(interp1(f,spectrum2(:,channel),[12:f:16]));
betahigh=mean(interp1(f,spectrum2(:,channel),[16:f:25]));

%creating an array
%store all iteration loop outputs in a matrix

m(kk,:)=[theta alpha betalow betahigh];

end

To get the total power across bins, sum the power in each bin. Also you need to compensate for your window loss if you want an accurate result.

For a rectangular window, the power in each DFT bin is given as:

$$P = \frac{Y[k]Y[k]^*}{N^2} = \bigg(\frac{|Y[k]|}{N}\bigg)^2$$

and $$P_{dB} = 20Log_{10}\bigg(\frac{|Y[k]|}{N}\bigg)$$

The total power over a range of such bins would be:

$$P_{k} = \frac{\Sigma_k(|Y[k]|^2)}{N^2}$$

Windowing Loss

When we apply a window, the coherent gain and non-coherent gain are affected differently. If the N samples are coherent (such as a frequency completely contained in a bin), the output grows at rate N (direct summation). If the N samples are non-coherent (spread across bins), the output grows at the rate of $\sqrt{N}$ (rms summation).

If we are only concerned with the signal power in each bin, we would use the coherent gain which is given as

$$G_c = \frac{\Sigma w[n]}{N}$$

where $w[n]$ are the window coefficients.   (For example the Blackman(30) window has a coherent gain of 0.406 or -7.83 dB):

Likewise non-coherent samples grow at $\sqrt{N}$. The non-coherent gain refers to the change in DFT output level of the non-coherent samples (such as white noise) due to the window function.

The normalized non-coherent gain is the rms value of the window weights given as:

$$G_{nc} = \sqrt{\frac{\Sigma w[n]^2}{N}}$$

The SNR change or processing gain of the window would be the ratio of the two, in dB this would be:

$$PG = 20log_{10}\bigg( \frac{G_c}{G_{nc}} \bigg)$$

If you are only interested in the signal power level and not signal to noise ratio, so you would simply compensate for the power computed by using $G_c$.

It may confusing to some how the power for a single tone can be effected differently by the window than the power in a signal that is spread evenly across multiple tones as described above. An intuitive explanation of this is that windowing the time domain function results in each bin in the FFT having an equivalent noise bandwidth (ENBW) that exceeds one bin. A rectangular window in comparison has an ENBW of exactly one bin, so when we sum the power of all the bins in the frequency domain, Parseval's theorem proves correct in that the total power from that summation is equal to the total power in the time domain signal. However if we window in time with something other than a rectangular window, the effect of the wider ENBW means that each bin in frequency is now reporting the power for that bin PLUS adjacent bins. So in this case when we sum the reported power over all bins, the total summation will over-estimate the actual noise power. In a case of a single tone, or tones that are spaced wider than the ENBW, each bin will accurately report the power for each tone when we have scaled the result by the recripocal of $\sum W_n$ where $W$ is the windowing funciton. (for the rectangular window this scaling results to $1/N$).


question related to something in karlin and taylor stochastic processes one text This question is essentially a question about something in Karlin and Taylor's Stochastic Processes One text in the spectral chapter. Since this is a DSP list, Karlin and Taylor may not be so popular so I will describe the setting which I think I've seen in other texts also.

They define the random variable $X_n =  A  \cos(\omega n) + B \sin(\omega  n)$.

So, $A$ and $B$ are zero mean normal random variables. $\omega$ is fixed at some frequency between $0$ and $\pi$. The process is on  the integers so the process start at some value of $-k$ and continues through to positive $k$. So $k$ can be any integer and $X_n$ represents a discrete stochastic process.

They then say that the $X_n$ process is stationary in the mean. Evidently this is so obvious that there's no need to explain it. I don't get it. Well, I sort of half get it.

A and B are  mean zero so, okay the mean will always be zero no matter what the values of the cos or sin term.  But since $n$ is increasing by 1 as the process proceeds, isn't the trajectory of $X_n$ changing at each $n$ because it multiplies omega ? So, because A and B are zero, stationary in the mean works. But, say A and B were both 1.0 say. Would the $X_n$ process still be stationary in the mean ? Since $X_n$ and $X_{n+1}$ will be at a different frequencies, $\omega  n$ and $\omega (n+1)$ respectively, doesn't that make the  mean different at times $n$ and $n + 1$ ? so it can't be a stationary process.

So, I guess my question is: For this process,  does stationary in mean only hold because $A$ and $B$ are zero in mean or would it work no matter what their means were ? If the latter is true, then I'm totally not understanding the process. Thanks.

The mean, i.e., the expected value of $X_n$ is

$$E\{X_n\}=E\{A\cos(\omega n)+B\sin(\omega n)\}\tag{1}$$

Since the expectation operator is linear, $(1)$ is equivalent to

$$E\{X_n\}=E\{A\cos(\omega n)\}+E\{B\sin(\omega n)\}\tag{2}$$

and since $\cos(\omega n)$ and $\sin(\omega n)$ are deterministic you finally get

$$E\{X_n\}=E\{A\}\cos(\omega n)+E\{B\}\sin(\omega n)\tag{3}$$

which equals zero if both $E\{A\}$ and $E\{B\}$ are zero. If at least one of the two random variables $A$ and $B$ has a non-zero mean, then $(3)$ will generally depend on $n$ (except for very specific choices of $\omega$), and, consequently, the mean of $X_n$ will also depend on $n$, i.e., $X_n$ cannot be stationary.

As a side note, a constant mean does of course generally not imply stationarity of $X_n$. For the given process, however, if $A$ and $B$ are identically distributed and $E\{A\}=E\{B\}=0$, and under the additional assumption that $E\{AB\}=0$ holds, $X_n$ is a wide-sense stationary (WSS) process.


How do I compute the gradient vector of pixels in an image? I'm trying to find the curvature of the features in an image and I was advised to calculate the gradient vector of pixels. So if the matrix below are the values from a grayscale image, how would I go about calculating the gradient vector for the pixel with the value '99'?

21 20 22 24 18 11 23
21 20 22 24 18 11 23
21 20 22 24 18 11 23
21 20 22 99 18 11 23
21 20 22 24 18 11 23
21 20 22 24 18 11 23
21 20 22 24 18 11 23

Apologies for asking such an open ended question, I've never done much maths and am not sure how to start tackling this.

So, if you're using Matlab, you can do:

X = your matrix

[gx,gy] = gradient(X); % first order gradient

[gxx,gxy] = gradient(gx); % second order gradient

[gxy,gyy] = gradient(gy); % second order gradient

To find the curvature of features I advise you to look into the eigenvalues and eigenvectors of a Hessian matrix. A hessian matrix is a square matrix of second-order derivatives of a scalar valued function. In this case the scalar field is the intensity, and the second order derivatives are $g_{xx}$, $g_{xy}$, $g_{yy}$.


Increase SNR by using only the first 200 samples I am making an experiment using sound source and microphones. I use metronome as the sound source, so I can consider it as impulse sound. I placed the microphones at different places about 1 meter above the floor. I use a recording system to get a sound file from each microphone. However I noticed that the acquired sound waveform from each microphone is different. For example let's say that one impulse sound is about 300K samples (48K sampling rate). The waveform of sound in microphone 1 and microphone 2 is similar only about the first 500 samples, mic 1 and mic 3 is similar only about the first 200 samples and so on. My guess that it is caused by the reverberation from the environment. Using Matlab function to find the delay is only good when I use only the first 200 samples.
Question:

Can I consider this as low SNR?
Can I use only the first 200 samples for each impulse sound as to improve the SNR?

Thank you.

If the SNR is limited by reflections in the room, and the impulse response of the channel from the source to a specific microphone is unknown, then the received signal can be used to determine the impulse response of the channel (together with a known transmitter signal that is spectrally rich in content) and with that equalize the received signal. However 200 samples would likely be far too insufficient for this purpose. In channel equalization we are basically solving for inverse convolution a system of overdetermined equations such that each output (sample) represents a sum of prior samples with unknown weights (which we solve for). These prior samples must extend over the significant duration of the channel impulse response, and then we need several samples to determine a least squared solution: each new sample provides us a new equation, and we would like to have more equations than unknowns as the accuracy of the solution will be improved the more equations we can use -- up to the limits of assumed stationarity in the channel statistics (the channel will change with time, so up to a point further inputs will contribute no further improvement, and then it will actually degrade the result).
200 samples at 48 KSps is only 4.16 ms. I believe typical channel impulse responses for a small room can be in the 100's of ms but I am less familiar with audio engineering to have confidence in that. But to equalize the channel to improve SNR by eliminating reflections/reverb I would recommend determining the impulse response time of the channel and then using at least 5 times more samples than that for the channel estimation with a known spectrally rich signal (white noise such as a maximum length sequence or frequency chirp). The equalization filter itself once the coefficients are determined would be the length of the channel response.
Once all channels are equalized for each microphone, assuming the remaining background noise is white and uncorrelated, further SNR improvement can be gained by coherently summing the microphone outputs: align their delay and sum each signal weighted by the SNR for that microphone for optimal ratio combining. If all microphones have the same SNR, this will provide a $\sqrt{N}$ SNR improvement where $N$ is the number of microphones.
If SNR is not limited by reflections in the room or other coherent sources, then any number of samples from multiple microphones can be added to get the $\sqrt{N}$ SNR gain assuming the coherent signals are aligned and the noise is uncorrelated on each microphone.


Calculating the power of a signal I'm trying to calculate the power of a signal and my tutor has given me this formula to do it

I've spent the past while building a program and now the foundations are there it's time to implement the maths side of it. The problem is I can't actually read it. Would someone be able to transcribe it and explain it to me?

As @jojek has already said in the comments, the formula reads
$$
P_\mathrm{x} = \dfrac{1}{T}\sum_{t=1}^{T}x^2(t)
$$
As $t$ usually denotes continous time I find this formulation a little bit odd. What your tutor probably actually meant is
$$
P_\mathrm{x} = \dfrac{1}{T}\int_{t=0}^{T}x^2(t) \mathrm{d}t
$$
It calculates the average power of the (continous-time) signal $x(t)$ of length $T$. In computer programs we're dealing with discrete-time (and discrete valued) signals, of course and the average power of the digital signal $x_n$ is given by
$$
\tilde P_\mathrm{x} = \dfrac{1}{N}\sum_{n=1}^{N}x_n^2,
$$
where $N$ is the length of $x_n$ in samples.


Sequence Length for a Linear Feedback Shift Register What is the total length for the sequence for a linear feedback shift register generated with a maximum length sequence, before repeating?

The duration of the linear feedback shift register (LFSR) in number of "chips" is $2^N-1$ where N is the number of states in the shift register (and order of the generator polynomial), and chips refers to each unique output that is generated in the sequence. This is assuming that you are using a maximum length sequence in the implementation of your shift register, which means the generator polynomial is both primitive and irreducible in GF(2).

Primitive: Repeated multiplying by x (representing a shift of one position in the shift register) will generate every possible state (except for the all 0 state). If the generator is not primitive, possible states (combinations of the register values) will be skipped and the sequence will be less than $2^N-1$ in length before repeating.

Irreducible: Cannot be factored, see example:

Reducible: $1+x+x^2+x^3 = (1+x^2)(1+x)$

Irreducible: $1+x+x^3$

Note the above polynomials are in GF(2) so the only elements are 0,1. 1+1=0, 0-1=1, etc

All primitive polynomials are irreducible but the converse is not true!

Primitive and irreducible polynomials are well tabulated. For a handy table of primitive and irreducible polynomials in GF(2), see Peterson's Table of Irreducible Polynomials

The graphic below shows the two common implementation approaches for the LFSR from the generator polynomial, in this case N=4, using generator polynomial $1+x^3+x^4$ which would generate $2^4-1=15$ chips before repeating:


Complex Least Squares FIR filter design using Levinson' s algorithm i'm trying to implement an arbitrary IIR filter design. I'm decomposing all the algorithm into smaller pieces that i can test.

I currently am on the part where it designs a FIR filter from an arbitrary response.

Matlab's main function is here :

function h = lslevin (N ,om, D, W)
% Complex Least Squares FIR filter design using Levinson' s algorithm
%
% h filter impulse response
% N filter length
% om frequency grid (0 <= om <= pi)
% D complex desired frequenccy response on the grid om
% W positive weighting function on the grid om
om = om(:);
D = D(:);
W = W(:);
L = length(om);
DR = real(D) ;
DI = imag(D);
 a = zeros(N,1);
 b = a;
% Set up vectors for quadratic objective function
% (avoid building matrices)
dvec = D;
evec = ones(L, 1);
e1 = exp(j*om);
for i=1:N,
    a(i) = W.'*real (evec);
    b(i) = W.'*real (dvec);
    evec = evec.*e1;
    dvec = dvec.*e1;
end
a=a/L;
b=b/L;
% Compute weighted l2 solution
h = levin(a, b) ;

and levin function :

function x = levin(a ,b)
% function x = levin(a ,b)
% solves system of complex linear equations toeplitz (a)*x = b
% using Levinson's algorithm
% a ... first row of positive definite Hermitian Toeplitz matrix
% b ... right hand side vector
%
% Author: Mathias C. Lang, Vienna University of Technology, AUSTRIA
% 9-97
a = a(:);
b = b(:);
n = length(a);
t = 1;
alpha = a(1);
x = b(1) / a(1);
for i = 1:n-1,
    k = -(a(i+1 :-1: 2)' * t)/alpha;
    t = [t;0] + k*flipud( [conj(t);0]);
    alpha = alpha*( 1 - abs(k)^2);
    k = (b(i+1) - a(i+1:-1: 2)' * x) /alpha;
    x = [x;0] + k*flipud(conj(t));
end

an example function can be done here (it should be a simple bandpass):

om1=pi*[linspace(0,0.23,230),linspace(0.3,0.5,200),linspace(0.57,1,430)];
D1=[zeros(1,230), exp(-j*om1(231:430) *2), zeros(1,430)]; % zeros are bandstop others bandpass
W1=[10* ones(1,230), ones(1,200), 10 *ones(1,430) ];%weight function
h = lslevin(61,om1, D1, W1) ;
wavwrite(h, 44100,  'e:\tmp2.wav'); % write the impulse response

but when i compute the Impulse Response, i have something like that :

which is far from the original design

can you help me pinpoint the errors in the code ?

thanks a lot

Jeff

EDIT1 : here is the final FFT of the IR from the answer given below:

From this term in D1:

exp(-j*om1(231:430) *2)

it looks like you're asking it for a filter with a group delay of 2 samples even though the number of taps is 61.  Try setting this to:

exp(-j*om1(231:430) * 61/2)

i.e. set a desired group delay of about half the filter length.


Implementing digital filter by padding with zeros I am making a 9th order digital bandpass filter with lower and upper corners of 200 kHz and 40 MHz respectively. I am using this filter to filter a 1D time domain signal which is 64k samples long sampled at a frequency of 150MHz.

I have done some digital filtering before in university, so I know what to expect, but it's been a while.

I have used this site: http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html
to generate the filter coefficients and the gain, and I can see how the code works, my question is this:

I start with the value of X[0], then to calculate the value of Y[0] I need values for Y[-1] ... Y[-18] and x[-1] ... X[-18].

I know these do not exist, so I think (from what I remember from university) that I pad with zeros, however, doing a bit or reading I heard it mentioned that padding with zeros changes the sampling frequency.

So how do I go about calculating this new sampling frequency (if indeed required..)?

Zero padding your input signal does not change your sampling frequency, it just changes (in this case probably very slightly) the time duration of your input signal.  An example would be that whether you filter one second's worth of data or two second's worth of data doesn't have any impact on the sampling frequency of the data.

One thing to keep in mind is that there will be a start up transient when zero-padding because your input signal instantaneously goes from zero to the value of the signal, essentially applying a step function on top of your signal of interest.  Generally speaking, you'll want to make sure you have some amount of data before a feature you're really interested in to make sure the start up effects have sufficient time to decay.  The website you mention gives the step response of the generated filter so you can check to see what effect it will have on your signal of interest.


Compute output given input, transfer function and initial conditions The problem statement is

Consider a causal LTI system whose transfer function $H(s)$ is given as $$H(s)=\frac{s+2}{(s+3)(s+4)}$$ Compute the output $y(t)$ for an input $x(t)=e^{-2t}u(t)$ when $y(0)=1$ and $y’(0)=0$.

The solution is given as $y(t)=5e^{-3t} − 4e^{-4t}$.
Here is my attempt:

Laplace transform left side and right side and consider the initial condition

right side is 0 , and left side is below

so i think answer is below but the prof's answer is not match with my answer.

Your solution looks correct. In any case, the given solution $y(t)=5e^{-3t}-4e^{-4t}$, $t>0$, must be wrong because it doesn't satisfy the initial condition on the derivative ($y'(0^+)=0$):
$$y'(t)=-15e^{-3t}+16e^{-4t},\qquad t>0\tag{1}$$
So we have $y'(0^+)=1$. Your solution satisfies both initial conditions.


How to extract noise when the noise on three data sets are near identical? I have a microphone array with three microphones placed in a small equilateral triangle. During recording, there is a lot of noise introduced by an ADC. However, I am not interested in reducing this noise before sampling - for learning, I want to remove this noise digitally as best as possible, given that the noise on each channel is near identical.

Given that their noise is near identical, I imagine it would be possible to correlate the signals in some way, even in the presence of a stronger signal with slight time delays. If possible, how can I remove this noise?

Given the three measurements $y1$, $y2$ and $y3$ coming from three unknown signals $x1$, $x2$ and $x3$, which we could assume additive to some unknown common noise $x4$, we can write:
$$
y1=x1+x4\\
y2=x2+x4\\
y3=x3+x4
$$
or equivalently:
$$
y=Ax\\
A=[1,0,0,1;0,1,0,1;0,0,1,1]
$$
Again, assuming a enough noisy signal, we can have an standard Recursive Least Squares algorithm for solving the standard LS projection recursively:
$$
\hat x=(A^TA)^{-1}A^Ty
$$
This way, and if and only if the additive noise and the noisy signal suficciencies are met, you will solve for your three signals, and your coupled noise.

Because of the former, you even can express this solution using the Kalman Filter


How to get Step Response from an Impulse Response in Simulink? Im sorry if it seems like a trivial question but I am confused right now. Im currently watching my professors lecture. And he basically has $\lambda(t)$ as Step response and $k(t)$ as Impulse response of an LTI system. The integral of the Impulse response gives us the Step response:
$$ \lambda(t) = \int^{t}_{0} k(t) dt $$
and
$$ k(t) = \lambda^{\prime}(t) $$
So to see it in action I decided to make a Simulink model of that and here it is (at least I think its good if not please tell me how to do it right):

As you can see those plots are not identical but I was expecting them to be (?)
Why is that?

Your system H(s) is continuous

$$ H(s) = \frac{s}{s+1} $$

Using a discrete impulse does not make a lot of sense. I suspect Simulink added a "zero-order hold" or ZOH block between your discrete impulse and your system $H(s)$. As such, your impulse response will be affected by the sampling period of your discrete impulse source. If you decrease the sampling period, the second curve shoud match the step response more closely, although they will never be exactly the same.

Edit :

Decrease the width of your impulse to say 0.1 second and increase the amplitude to 10 so that the area under the curve stays 1.


Image Contrast Enhancement based Sub-histogram Equalization Technique without Over-equalization Noise I'm a computer science student and I have an assignment requiring me to implement the algorithm proposed by this paper but I'm hitting a roadblock trying to understand the implementation of one of the methods proposed in the paper. I'm not sure if I understand Section IIIB well.
This is what I understand so far.

In order to separate the histograms into subregions, we must first obtain the histogram of the intensity levels. Once achieved, we calculate the density of the histogram (equation 7).
Furthermore we have to sort the density function values in descending order. How is this implemented? I assume its only in the following paragraph in the paper that goes on to describe this but I'm having a hard time understanding what the authors are trying to convey.

Once we find the groups of high density areas, we extend a group area surrounding these values. The paper talks about momentary changes. I'm not sure what this means, anyone have an idea?

The page is located at
http://www.waset.org/journals/ijeee/v3/v3-6-53.pdf

The paper is not really well written, so it's not surprising that you're having a hard time. These are a few hints to guide you:

the authors try to do some mode detection on the image brightness histogram in order to segment it between common values in the image (= high density in their terminology) and rare values in the image (= low density);
in order to segment the histogram, they fix a window width (i.e., a number of bin of the histogram), apply this window over all possible positions in the histogram and sum the populations of the bins inside the window ($Eq. (7)$). This gives the density of a histogram bin (or equivalently of a brightness level);
their segmentation process consists then in going from the peaks of the histogram and merging neighbouring bins into it as long as the difference between the border bin and its neighbour (= the momentary change) is less than the average of the momentary changes in the region.


Calculate ifft using only REAL forward FFT There are already questions on this topic here and here. I am trying to calculate the ifft using a DSP that only has a library function for a forward REAL (assumes conjugate symmetry) fft. For a given a signal of length N, the output of the FFT process is length N interleaved real/complex values, where the real part of the nyquist bin is put into the imaginary part of the DC component (always zero for a real signal - so this keeps a constant memory footprint). I'm having trouble fathoming how to calculate the inverse using such a function.

EDIT: This answer is basically a more detailed version of Hilmar's answer, which I noticed too late. Below there are more details about the derivation and the rearrangement of the outputs of the real FFT.

I assume you use the definitions of the FFT and IFFT as commonly used in DSP:

$$\textrm{FFT: }\quad X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

$$\textrm{IFFT: }\quad x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{2}$$

For a real FFT we assume that $x[n]$ is real-valued. We also assume that the FFT length $N$ is even (otherwise you don't get an FFT bin at Nyquist).

From $(1)$ and $(2)$ we see that

$$\text{IFFT}\big\{X[k]\big\}=\frac{1}{N}\left(\textrm{FFT}\big\{X^*[k]\big\}\right)^*\tag{3}$$

Since your FFT routine only accepts real-valued inputs, we need to write $X[k]$ in terms of its real and imaginary parts:

$$X[k]=X_R[k]+jX_I[k]\tag{4}$$

Furthermore, if we assume that the inverse FFT of $X[k]$ is real-valued, we can rewrite $(3)$ as

$$\begin{align}\text{IFFT}\big\{X[k]\big\}&=\frac{1}{N}\textrm{Re}\Big(\textrm{FFT}\big\{X_R[k]\big\}-j\,\textrm{FFT}\big\{X_I[k]\big\}\Big)\\&=\frac{1}{N}\Big[\textrm{Re}\big(\textrm{FFT}\big\{X_R[k]\big\}\big)+\textrm{Im}\big(\textrm{FFT}\big\{X_I[k]\big\}\big)\Big]\tag{5}\end{align}$$

A little extra inconvenience lies in the fact that a real FFT routine only computes $N/2+1$ values whereas Eq. $(5)$ assumes FFT results of length $N$. This Matlab/Octave script shows how to correctly arrange the results of the individual FFTs:

N = 64;     % must be even!
N2 = N/2;

x = randn(N,1);

X = realfft(x);
XR = real(X); XR = [XR;XR(N2:-1:2)];
XI = imag(X); XI = [XI;-XI(N2:-1:2)];

x2a = real(realfft(XR)) / N; x2b = imag(realfft(XI)) / N;
x2 = [x2a + x2b; x2a(N2:-1:2) - x2b(N2:-1:2)];

In the above script we assume that there's a function realfft.m defined by

function X = realfft(x)

N = length(x);
if rem(N,2), error('N must be even.'); end

X = fft(x,N);
X = X(1:N/2+1);

end


Sawtooth wave Fourier coefficients I have to calculate the Fourier coefficients of this signal. I found that signal equation is $$ y = \frac {A(2t-T)}{T} $$

To find Fourier coefficients I wrote

$$ x_k = \frac{2A}{T} \int_{0}^{T/2} \frac{2t-T}{T} e^{-i2 \pi k f_0’ t } $$

In this case $$ f_0’ = f_0 $$ because the period is T_0

I calculate the integration by parts of the first integral and I obtained $$ - \frac{-e^{-i \pi k }( i \pi k ) +1 - e^{-i \pi k }}{4 (\pi k f_0 )^ 2} $$

( i control This with wolfram and should be correct, I hope !)

After solving the integral by parts I obtained that

$$ x_k = \frac{A( e^{-i \pi k} + 1 - \pi i k )} { (\pi k )^2} $$ but the result should be $$ \frac{iA}{k \pi} $$

I can also post all my resolution. Can someone help me ? Thank you

This is one of the best problems to demonstrate Fourier Series properties, and specifically the time derivative property: $$\frac{d}{dt}x(t) \overset{FS}\longleftrightarrow j2\pi kf_0 X_k$$

Instead of computing the integral $$X_k=\frac{1}{T_0}\int_{T_0}x(t)e^{-j2\pi kf_0t}dt$$ which is time consuming, you can take the first derivative of your signal and compute the Fourier coefficients of the derivative. At discontinuities, you will have Delta functions, and everywhere else you will have a constant equal to the slope of your $y$ signal, which is $2A/T_0$. In other words, inside the time interval $[0,T_0)$, your derivative can be expressed as $$\frac{d}{dt}x_{T_0}(t) = \Bigg\{\begin{array}{ll}\frac{2A}{T_0}, & t \in (0,T_0) \\ -2A, & t=0\end{array} = \frac{2A}{T_0}\mathrm{rect}\Big(\frac{t-\frac{T_0}{2}}{T_0}\Big)-2A\delta(t)$$
I will not use the last equality for integration, but you now see it is much easier to integrate $\frac{d}{dt}x_{T_0}(t)$. Let's say $z(t)=\frac{d}{dt}x(t)$, so its Fourier coefficients would be $$\begin{align}
Z_k &= \frac{1}{T_0}\int_0^{T_0}\frac{2A}{T_0}e^{-j2\pi kf_0t}dt -\frac{2A}{T_0}\int_{0}^{T_0}\delta(t)e^{-j2\pi kf_0t}dt \\
&= \frac{2A}{T_0^2}\frac{1}{(-j2\pi kf_0)}e^{-j2\pi kf_0t}\Big|_{0}^{T_0} - \frac{2A}{T_0}e^{-j2\pi kf_0t}\Big|_{t=0} \\
&= \frac{2A}{-j2\pi kf_0T_0^2}(e^{-j2\pi kf_0T_0}-1)-\frac{2A}{T_0} \\
&= 0 - \frac{2A}{T_0} = -\frac{2A}{T_0}
\end{align}$$
Can you apply the time derivative property to find your $X_k$ coefficients?

P.S: If you are not aware of the Delta function and its properties, or even the properties of Fourier Series, then just skip my answer. :)


Detecting and removing interferences from a signal I am using MATLAB in order to denoise and remove interferences on a signal.
I used wdenoise to denoise my signal which works by setting a threshold (for example SURE) for each scale and set all coefficients below this threshold to zero (these coefficients represents noise). It works pretty well.
I also wanted to detect interferences, which are in my case like spikes (most of the time with high amplitude compared to the useful signal), and remove them. I try to use wavelet decomposition (wavedec) to detect theses spikes and set a threshold to identify them for each scale. If a coefficient is greater than this threshold, I consider that this is an interference coefficient so I set it to zero. Let's consider a decomposition of 8 level. Because of downsampling, make a threshold on D{8} for example, give me a signal with only 10 coefficients so it's very complicated to see what coefficient is an interference coefficient.
With this technique, I don't have great results. Any idea to reconstruct a signal with spikes? Maybe other techniques ? Maybe used swt or modwt instead of wavedec ?

I'd recommend using a median filter to smooth your signal. A median filter will get rid of the outliers, or the spikes in your signal. The length of the median filter will have to be determined by how frequently the spikes appear in your signal. The median filter is available in Python as numpy.median.
If you want to retrieve the spikes (the interference signal), you can simply subtract the filtered signal from the noisy signal to get the interference.


Difference between repeated poles and distinct poles? An important concept in dsp is marginal stability where we often see the term" repeated roots " or "repeated poles "?
What are they?

Does the term repeated means that two or more poles occur at exactly same coordinates?

Please kindly give answer along with examples for both s domain and z domain?

Repeated poles simply means there are more than one pole at the same location. If a pole is not repeated then it is a Distinct Pole. Consider the simple case of a cascade of two integrators (in s) which is similar to the cascade of two accumulators (in z).

The Laplace transform of an integrator
Is $\frac{1}{s}$, which is one pole at the origin. Two integrators in cascade would be $\frac{1}{s^2}$ so has repeated poles at the origin. The same with two accumulators in cascade would be $\frac{1}{(1-z^{-1})^2}$, which has repeated poles at z=1.

If there were two poles at z=1 and a third pole at z= 0.5, then such a system would have two distinct poles.


Resampling time series to regular array, then downsampling (Butterworth) Long time reader, first time poster.  I have a few very simple questions that are troubling me and I am hoping that one of you guys can help me out.

Setup & Aim: I have a time series that I want to downsample, and I simply want to run a lowpass filter on it before doing so to avoid aliasing. I am using Python (SciPy) but it looks like MATLAB behaves similarly, neither are really relevant for these questions.

My original time series is sampled at $0.5\textrm{ ms}$ ($2000\textrm{ Hz}, f_{\rm Nyquist}=1000\textrm{ Hz}$) and I want to resample to $2\textrm{ ms}$ ($250\textrm{ Hz}, f_{\rm Nyquist}=250\textrm{ Hz}$), so I must apply an anti-alias filter that cuts off any frequencies $> 250\textrm{ Hz}$, and then downsample. So far, so good.

In Python, it looks like a Butterworth Filter is the way to go, which requires a normalised frequency $\omega_n$.  My understanding is that in my case $\omega_n = 250\textrm{ Hz}/1000\textrm{ Hz} = 0.25$.

Now, what I don't understand and I cannot find any information on, is as follows:

What if my original time series ($f_s=2000\textrm{ Hz}$) had been upsampled from $1\textrm{ ms}$ ($f_s=1000\textrm{ Hz}, f_{\rm Nyquist}=500\textrm{ Hz}$)? There is no extra information between $500\textrm{ Hz}$ and $1000\textrm{ Hz}$ but I don't necessarily know that and I apply a Butterworth Filter with $omega_n = 0.25$ (instead of $\omega_n = 0.5$ for $1\textrm{ ms}$ sampling).  Is it an issue? Am I misunderstanding how a Butterworth filter works?
My second question is something like "Why is this the preferred implementation of a low pass filter?" I am sure there are good reasons but I have used software in the past to just high cut filter my data knowing my new $f_{\rm Nyquist}$, so in my case.  So in my case I would use something like $0-0-200-250\textrm{ Hz}$.  Again, what am I missing? I know that maximum frequency that I want to keep.

Finally, one of the roots of my problem is that I sometimes have irregularly sampled data.  I can run an interpolation to a regular time array but when doing this I tend to oversample, to avoid losing signal  (This is where my first question comes in).  Am I wasting my time? Should I just resample to the smallest time interval in my data?

I am not sure why you concluded "A Butterworth Filter" is the way to go, and looks like you have come to the same question- what originally led you to that?

A Butterworth Filter implementation stems from copying analog filter techniques and from what I have learned with DSP is that you have the wonderful opportunity to NOT copy the analog. This is because analog design techniques have been derived from the mathematical limits of what standard R L and C components can give you, while in DSP you can more likely address the problem from an optimum mathematical perspective within the greater space of delays, adds and multiplies (to note and hopefully not confuse; delays, adds and multiplies are NOT strictly limited to digital implementations; for example SAW filters). This may be a good case in point if we look closer at why you would choose a Butterworth filter versus what optimized algorithms such as least squares (FIRLS in Matlab) can give you. The least squares algorithm will result in a linear phase FIR solution, with a stop band roll-off which is ideal for multi-rate resampling applications. fred harris has made a good case (I cannot find the specific reference) that for resampling applications, and FIR filter will always out-perform an IIR filter.

Further what is really great about FIRLS for multirate applications is that you can easily implement a multi-band filter, concentrating the rejection bands to be just at the alias locations, as in the example shown in the figure below showing the same filter design for either a decimation (DFIR) or interpolation (IFIR) design, in this case using a factor of 4:

To answer your question about upsampling: If the signal has been upsampled, this is typically done by inserting zeros (which creates all the aliases in frequency) followed by a low pass filter which will eliminate the aliases and in the process grow the zeros to their optimum interpolated value, so yes in this case, as long as the intepolation filter is present, you will not need an additional low pass filter.

I show this in the following figure showing that if you do the interpolation first in a rational rate conversion (both interpolate and decimate) you can combine the IFIR and DFIR designs as one filter.

Also for resampling given its simplicity, consider using a CIC structure. If you are concerned about passband roll-off- the inverse Sinc compensator is quite simple (3 taps). Both are detailed in the following links:

Decimation concepts: Fast Integer 8 Hz 2nd Order LP for Microcontroller

Upsampling Spectrum: Upsampling and comb function

Inverse Sinc Filter: how to make CIC compensation filter

More about CIC: CIC Cascaded Integrator-Comb spectrum


How to remove the noise without destroying the main edge? I am working on a project. I want to segment a logo from a car and the picture's size is 3648*2432. The logo is selected by the red rectangle.

I extracted this area and turn into gray level. The area's size is 249*173. Then, I use histogram equalize to enhance the edge. Finally, I use adaptive canny to find edge.

However, the result is not perfect.
So, I used median filter to smooth the left picture. The result is as following.

I remove the small parts.

The result has been improved, but it is also not good. Besides, I have tried  bilateral filter, mean filter and Gaussian filter, but the result is not good, too. I know that after passing histogram equalize, the noise can be enhance. But, in some case, low contrast, I need to adapt this step to enhance my picture. Is any one who can provide me any different comment? Thank you very much.

Use bilateral filter or anisotropic diffusion first.

The effect of anisotropic diffusion is as the following:

. The MATLAB code can be found here.

Here is its effect on your image:

Finally, non-local means is a also a good way to get rid of the noise. You might also want to take a look into that. I warn you though, it is slow.


SIFT About Difference-of-Gaussian function extrema? How to get formula (2) by formula(1)?

You have to take the derivative with respect to the vector $x$ and set it equal to zero. For a constant matrix $A$, the derivative of $A^Tx$ is $A$, and the derivative of $\frac12 x^TA^Tx=Ax$. So taking the derivative of $(1)$ gives

$$\frac{\partial D}{\partial x}+\frac{\partial^2D}{\partial x^2}x\tag{1}$$

Setting $(1)$ equal to zero results in $(2)$.


Why does an fftshift in the time domain give a different result compared to the frequency domain? I am trying to speed up some of my code and for physical reasons I need to do an FFT shift to get the FFT data in the standard DFT fourier representation. In Matlab this is performed in the frequency domain, with a data swap.

I had believed that I found a numerical trick, by performing the shift in the time domain, before transformation I was able to save a reasonable amount of computation time.

The problem is that for some reason the results are slightly different. I am curious as to why.

The Matlab code is as follows for a 1K signal

%1khz
f=1000;
t=0:.0001:5;
x=sin(2*pi*f*t);
%Wave
subplot(4,1,1);
plot(x)
%FFT
subplot(4,1,2)
ffted=fft(x);
plot(real(ffted));
%FFT Shifted
subplot(4,1,3)
shifted=fftshift(ffted);
plot(real(shifted));
%FFT Shifted by preshifting
subplot(4,1,4)
s=zeros(1,length(x));%In languages like C++ this is computationally inexpesnive
for i=1:length(x)
    if (bitget(i,1))%odd evens, works for real signal
        s(i)=-1;
    else
        s(i)=1;
    end
end
chirpshifted=fft(s.*x);
plot(real(chirpshifted));
%Difference
A=sum(real(chirpshifted));
B=sum(real(shifted));
A
B

There are a couple issues with your code. First, let's talk about why you can do what you're trying to do.

What you're trying to take advantage of is the duality between a multiplication in the time domain and convolution in the frequency domain. The fftshift() function simply swaps the two halves of a vector (which is assumed to be the result of a call to fft()). This operation is equivalent to a circular rotation of half of the vector length in the frequency domain.

One can accomplish the same operation by multiplying by a complex exponential function in the time domain before the FFT, as you noted. To get an $M$-sample shift, you need to multiply by a complex exponential that looks like:

$$
s[n] = e^{\frac{j2\pi nM}{N} }
 $$

For the case where $M = \frac{N}{2}$, this simplifies to:

$$
s[n] = e^{j\pi n} = [1, -1, 1, -1, \ldots]
 $$

This illustrates problem #1 with your script: it multiplies by the inverse of the above: $-s[n]$. You need to invert the sense of the if statement in the for loop to ensure that the vector s has the appropriate value.

Problem #2 stems from another simple MATLAB mistake: all of the vectors in your simulation are length 50001. So, the time-domain rotation that you're doing by multiplying by the $s[n]$ illustrated above actually results in an effective frequency domain shift of 25000.5 bins (instead of the 25000-bin shift implemented by fftshift()). The mismatch in results is due to this extra half-bin shift.

You can fix this problem in two ways:

Change the time-domain multiplication to multiply by an actual complex exponential vector at the appropriate frequency (i.e. set $M=25000$ and $N=50001$ in the first equation). The resulting signal at the input to the FFT will be complex.
Change your simulation to use even-length vectors (i.e. chop off the last element of the t vector that you initialize at the beginning).

I tweaked your simulation to fix the generation of s and truncate the vector length to 50000. It is as follows:

%1khz
f=1000;
t=0:.0001:5;
t=t(1:length(t)-1);
x=sin(2*pi*f*t);
%Wave
subplot(4,1,1);
plot(x)
%FFT
subplot(4,1,2)
ffted=fft(x);
plot(real(ffted));
%FFT Shifted
subplot(4,1,3)
shifted=fftshift(ffted);
plot(real(shifted));
%FFT Shifted by preshifting
subplot(4,1,4)
s=zeros(1,length(x));%In languages like C++ this is computationally inexpesnive
for i=1:length(x)
    if (~bitget(i,1))%odd evens, works for real signal
        s(i)=-1;
    else
        s(i)=1;
    end
end
chirpshifted=fft(s.*x);
plot(real(chirpshifted));

I then compare the results of the two methods:

>> diff = abs(chirpshifted-shifted);
>> sum(diff)
ans =

   3.3406e-21

As you can see, the results agree as perfectly as you would expect from floating-point math.


Will signal get distorted if channel bandwidth is lesser than the signal bandwidth? This statement is given in the book "digital communicatinon by Taub schilling "

Now one point I am unable to understand that if the channel bandwidth is lesser than the signal bandwidth,how it is possible that the signal will not get distorted? can someone explain this to me in simple word Please?

The key here is the title of the book you are reading: "digital communicatinon by Taub schilling " as you say

What the book is talking about is digital communications and not analog communications in which the fidelity of the transmitted signal is of great importance; we want the output signal to be as close a replica of the input signal as is possible. In digital communications, it does not matter diddly-squat if the transmitted waveform is distorted as long as from the received waveform one can determine which of the two (more generally, few) possible transmitted waveforms $s_0(t)$ and $s_1(t)$ caused this received signal to occur. If the receiver can make this determination with high accuracy, then we have managed to transmit one bit (or more generally $\log_2$(few bits)) with high accuracy by choosing $s_0(t)$ or $s_1(t)$ according as we wish to send a $0$ or a $1$. Even if the received signal
$r(t)$ is badly distorted from what was sent, the receiver can still determine whether the transmitter sent $s_0(t)$ or $s_1(t)$: fidelity of transmission is unimportant.

With the above in mind, note that the Taub and Schilling description is dreadfully ambiguous, and in many spots downright incorrect. An equalizer (a concept from digital communications) cannot perfectly undo the distortion caused by forcing a 1 kHz analog signal through a low-pass filter with 3 dB bandwidth $1$ Hz. An inverse (analog) filter is trivial to write down --- it is a filter whose transfer function is $\frac{1}{H(f)}$ where $H(f)$ is the channel response --- but in most cases that inverse filter is unrealizable. The truth of the matter is that any filtering by a channel causes the transmitted signal to be spread out in time and reversing this time spreading and bringing the signal back to its pristine shape is in almost all cases impossible: once the genie is out of the lamp, he cannot be stuffed back in, ditto Pandora's box for those preferring a more Western analogy. That last sentence "The signal is then recoverable precisely as transmitted" is sheer nonsense.


Is there an intuitive explanation for what complex basebanding a signal is? I have an algorithm which takes an input signal, $x[t]$, and (among other things) complex basebands the data and outputs the phase of the complex numbers.
I'm new to the field and the Matlab function which produces the basebanded signal is a kind of black box to me. I've tried reading the code, but it's pretty impenetrable and not commented. I also tried to input different shapes of signals, but couldn't learn much from it. My questions are:

What exactly is the complex baseband of a signal? I sort of understand what "baseband" means, but where does the complex phase come from?

Related, does the phase have any relation to the signal? If I input a sinusoid, what should I expect the phase output to be?

What exactly is the complex baseband of a signal? I sort of understand what "baseband" means, but where does the complex phase come from?

A baseband signal is one where the center frequency is $\approx 0$ Hertz. There may be some small offset, but the idea is that it is no where near a passband signal whose center frequency is in the millions or billions of Hertz.
For a real signal, the frequency spectrum is symmetric. Due to the up-conversion process (converting baseband to passband so that the center frequency is no longer $0$ Hertz), we can put data on two orthogonal carriers, $\text{cos}$ and $\text{sin}$, in the passband signal:
\begin{align}
x(t)&=\text{Re}\bigg(s(t)e^{j2\pi f_ct} \bigg) \\
&= x_I(t)\text{cos}\big(2\pi f_ct\big) - x_Q(t)\text{sin}\big(2\pi f_ct \big)
\end{align}
where $s(t)=x_I(t) + jx_Q(t)$, is the complex baseband signal. For a nice explanation of what is probably happening in your MATLAB function on how to convert passband to baseband, check out this previous answer: https://dsp.stackexchange.com/a/43281/31316.

Related, does the phase have any relation to the signal? If I input a sinusoid, what should I expect the phase output to be?

There is a relation...it is the phase, that is all. The phase is helpful though and can help to pack more bits of information into a single symbol.


Why not use a complex number in the exponent in the z-transform? I am trying to comprehend how the z-transform has come to be similar but different from its continuous counterpart, namely the Laplace transform.

It seems to me that the most parallel and intuitive approach to the deduction of a generalized form of the DTFT would be to substitute $j\omega$ for $z = \sigma + j\omega$, i. e.:
$$
X(z) = \sum_{n = -\infty}^{+\infty} x[n]z^{-n} = \sum_{n = -\infty}^{+\infty} x[n]e^{-(\sigma + j\omega)n}
$$

Instead, we find normally that $z$ is a polar-form complex number which is the product of its magnitude and $e^{j\omega}$. I can guess that probably such representation leads to a rectangular representation of the z-transform, eventually rendering the analysis of poles and zeros impractical due to the periodic nature of the DTFT, but I am not sure. Could anyone enlighten me?

$z^{-1}$ is the unit delay operator.

$z$ is the unit advance operator.

$s$ is the derivative operator.

$1/s$ is the integration operator.

the utility of the Z and Laplace transforms is to solve (analyze, predict) constant coefficient differential and difference equations algebraically.


why do directional derivatives work better with a gaussian filter? I've been working on image segmentation and computing the directional derivatives of a grey-scaled image, with the objective of detecting contours and edges. I have realised that if I apply a gaussian filter before computing the directional derivatives, the edges get enhanced better. So my question is - why does that happen?

Here's the original image:

This is what the first order derivatives (dx and dy) look like (without gaussian filter):

here's the derivatives after applying gaussian filtering to the image, with sigma = 6:

*******Python code to compute the derivatives

with rio.open("vein_template.tif", 'r') as ds:
    RGB_arr = ds.read(masked=True)  # read all raster values
rgb_1 = np.rot90(np.transpose(RGB_arr, [1, 2, 0]))
#the rgb image contains 4 channels, so I'm changing their order, and making the image upright

%matplotlib
plt.imshow(rgb_1)

# Transform to grey-level img
gsi = 0.2989 * rgb_1[:, :, 0] + 0.5870 * rgb_1[:, :, 1] + 0.1140 * rgb_1[:, :, 2]
plt.imshow(gsi, cmap = "gist_gray")

#compute derivatives
dx = np.diff(gsi, axis = 0)
dy = np.diff(gsi, axis = 1)

f = plt.figure()
f.add_subplot(1, 2, 1)
plt.imshow(dx, cmap = "Greys")
plt.title("dx")
f.add_subplot(1, 2, 2)
plt.imshow(dy, cmap = "Greys")
plt.title("dy")

Well, look at your original picture: it's constant for all points but the edges, which means your derivative is zero for all points but these edges.

By applying a "rounding, smoothing" filter to it, you "smear" the edges enough to make the derivative be non-zero for multiple pixels, in every direction.


Estimate $s$ from $y=s+n$ Now I receive $y=s+n$, $s$ is the signal and $n$ is the white guassion noise. Everytime $s$ will be change irregularly, have some methods to get a estimation of $s$ namely $\hat s$?

Namely $s$ is random.
And I don't know $s$ how to change everytime,so the traditional kalman filter algorithm may can't be used in this circumstance.

I think you are seeking for an algorithm but this may help. If you don't know how $ s $ is distributed but you know its moments of first and second order, you can try with the linear estimator of the linear Minimum-mean Square Error (lmmse), with the following form,
$$
\hat{s}(y) = ay + b
$$
with
$$
a = \frac{ \mathrm{C}_{xy} }{ \sigma_{y}^{2} } \\
b = m_{s} - am_{y}
$$
where $ m $ denotes mean, $ \sigma $ variance and $ C $ covariance. It would help to know the moments of the noise, I mean, for example if is zero-mean and normalized to unity.