Wavelength filter with OpenCV I am new to OpenCV and DSP. I want to code a filter which act like a real optical filter. When you filter the image, the output looks like shotted with a $550nm$ optical filter lens. Like a Bandpass filter for image. I wanna filter all lights below $550nm $ wavelength.

Is it possible with OpenCV? How can Ii do that? Filter2D or Gabor Filter?

Sorry for my English and sorry if the question is irrevelent.

Thanks

In general this is not possible, since you do not have full spectral information for each pixel, but only coordinates in some low-dimensional color space (i.e. RGB).

If your image data is in RGB, then the best you can do is to simply take the green color channel (which will be most responsive to 550 nm) and throw away the red and blue.

Depending on the actual response of your RGB channels, you might be able to better estimate the 550nm contribution via a linear combination of red and green (see Spectral Sensitivity).


What is the value range of log magnitude spectrograms? I'm quite new to this field and have a probably rather basic question:
Is there an upper bound for the values in a log-magnitude spectrogram generated like this:

log1p(abs(stft(signal) ** 2)

Obviously, the lower bound is 0 but is there also an upper bound (based on window size/hop length of the STFT)?

Let us consider a standard discrete "spectrogram" or short-term Fourier transform with a window $w$ and  frames of length $N$. Let $\|x\|_r = (\sum_{n=0}^{N-1}|x_n|^r)^{1/r}$ denote the $r$-norm of $x$ ($r\ge1)$. Then for each frame of $x$, the first windowed Fourier transform is:
$$F(\omega)=\sum_{n=0}^{N-1}x_n w_n e^{-j\omega n}\,.$$
Hence, by the Rogers-Hölder's inequalities, for all $\omega$, and any couple  $(p,q)$ such that $\frac{1}{p}+\frac{1}{q}=1$:
$$|F(\omega) |\le \|w\|_p\|x\|_q\,.$$
In particular, for classical positive unit windows, $\sum_{n=0}^{N-1} w_n =1 = \|w\|_1$. Choose $p=1$ and $q=\infty$, and you get (over all frames):
$$|F(\omega) | \le \max_n | x_n|\,.$$
Not the tighest upper bound, but simple enough. And  honestly, using Rogers-Hölder was seen a bit of an overkill here, since simply:
$$\left|\sum_{n=0}^{N-1}x_n w_n e^{-j\omega n} \right| \le \sum_{n=0}^{N-1}\left |x_n w_n e^{-j\omega n} \right| \le \max_n | x_n|\sum_{n=0}^{N-1}\left | w_n  \right|\,, $$
but for "unit-energy" windows, you will get other estimates with $p=q=2$.
You can now plug it inside your "log of one plus" formula.
Similarly, you could find a tighter lower bound.


Detecting a low frequency cue in audio recording I am hoping to find a linux program or library for a language like perl or python (open to any solution though) to detect low frequency cues in an audio recording.

Ideally I would end up with a printout of timestamps in HH:MM:SS format that correspond to where the cues occur.

Here is visual example of the audio cue:

The cue appears to be unique waveform and there is no other sounds at this frequency through out the recording. The frequency is so low that it is in audible. Based on this I assume it must be possible to detect the cues.

You can solve this problem with a simple time based low pass filter, no need for an FFT:

Calculate the minimum integration time you want by dividing your sampling rate by the frequency you want to detect. For a 10Hz at 44100Hz, it would be 4410 samples, so let's take 4500 samples.
Then apply a low pass filter to this group of samples (choosing the right frequency cut), you will isolate the frequencies you want. I do not really know how to use them, but it should look something like this:

$$y(n) = x{(n)}+x(n-1) + ...$$

Calculate the relative power of your filtrated group of samples in dBFS (Full Scale) based on this formula:

$$p_{RMS} = \sqrt{\frac{ x_1^2 + x_2^2 + \ldots}n } $$
$$dbFS = 20\log_{10}\frac{ p_{RMS}}{p_{max}}$$
'n' behind the number of samples and 'pmax' the maximum value of a sample.

If the level you get is superior to a specific threshold of your choice (not to low, because your audio certainly have a tiny low frequency component), then there is a cue, otherwise there isn't.

As for the timestamp aspect, you can increase a variable each time you get into this low frequency detection loop, and calculate the time passed using the sampling rate. Then write an array containing the input/output time of each successful detection.

Not sure if the minimum integration time is really necessary though.


Realization of a filter based on its transfer function How can we check whether the filter is realizable given its transfer function and What are the parameters the realization depends on?
Here is an example:

Show that a filter with transfer function
$H(ω)=\frac{2(10^5)}{ω^2 + 10^{10}} e^{-jωt_0 }$
is unrealizable. Can this filter be made approximately realizable by
choosing a sufficiently large t0.

--Modern digital and analog communication systems (3rd edition) - B.P. Lathi
Note: The question is of general nature. The example quoted above serves only as to make it clear, what is being asked. No one has to provide a solution for this problem.

A transfer function is called realizable if it can be implemented by a causal and stable system. The given frequency response is continuous and doesn't have any impulses, so the corresponding system is stable.

The transfer function (as a function of $s$) is given by

$$H(s)=\frac{a}{b-s^2}e^{-st_0}\tag{1}$$

The given frequency response is obtained form $(1)$ by substituting $s=j\omega$. Since $b>0$, $H(s)$ has two real-valued poles at $\pm\sqrt{b}$, and, consequently, the corresponding impulse response is two-sided (i.e., non-causal).

Partial fraction expansion of $(1)$ gives

$$H(s)=\frac{a}{2\sqrt{b}}\left[\frac{1}{\sqrt{b}-s}+\frac{1}{\sqrt{b}+s}\right]e^{-st_0},\qquad -\sqrt{b}<\textrm{Re}\{s\}<\sqrt{b}\tag{2}$$

From $(2)$ it is straightforward to obtain the impulse response

$$\begin{align}h(t)=\mathcal{L}^{-1}\{H(s)\}&=\frac{a}{2\sqrt{b}}\left[e^{\sqrt{b}(t-t_0)}u(-(t-t_0))+e^{-\sqrt{b}(t-t_0)}u(t-t_0)\right]\\&=\frac{a}{2\sqrt{b}}e^{-\sqrt{b}|t-t_0|}\tag{3}\end{align}$$

No matter how large you choose $t_0$, the impulse response will never be causal. However, since the system is stable, the impulse response decays for $|t|\to\infty$, so you can approximate the actual impulse response well by sufficiently shifting it to the right (i.e., by choosing some large $t_0$), and truncating the non-causal part.

Also take a look at this related question and its answer.


Why do we do oversampling in raised cosine? I would like to do know why do we need to do oversampling in a raise cosine filter. There is already a question like this but the answer is very confusing and not sufficient.

I need to know what this oversampling factor does?
The higher is better or the lower?

A (baseband) pulse amplitude modulated (PAM) signal is given by

$$s(t)=\sum_{k}A_kp(t-kT)\tag{1}$$

where $A_k$ are the data symbols, $p(t)$ is the transmit pulse, and $1/T$ is the symbol rate. Note that $(1)$ describes a continuous-time signal. If you want to simulate a digital communication system, you need to use sampled versions of continuous-time functions, such as the transmit pulse $p(t)$. For this purpose, you need to define a sampling rate greater than the symbol rate. The ratio between the sampling frequency and the symbol rate is the oversampling factor.

Assume you choose a sampling rate an integer (oversampling) factor $m$ greater than the symbol rate:

$$\frac{1}{T_s}=\frac{m}{T}\tag{2}$$

From $(1)$, the sampled signal $s(nT_s)$ is then

$$\begin{align}s(nT_s)&=\sum_{k}A_kp(nT_s-kT)\\&=\sum_{k}A_kp\left(\left(n-k\frac{T}{T_s}\right)T_s\right)\\&=\sum_{k}A_kp((n-km)T_s)\end{align}\tag{3}$$

Or, in notation for discrete-time signals

$$s[n]=\sum_{k}A_kp[n-km]\tag{4}$$


formula derivation for periodic signal power Studying DSP on my own. Taking an introductory class. Was given a formula for signal power

$$ P = \lim_{M\to\infty} {\frac{1}{2M+1} \sum_{n=-M}^M} \lvert x[n] \rvert^2 $$

I do understand why the formula looks the way it does.

But then we were given a formula for a power of a periodic signal, where $N$ is the period,

$$ P = {\frac{1}{N} \sum_{n=0}^{N-1}} \lvert x[n] \rvert^2 $$

I can't figure out how this formula is derived from the original one. How do I do it?

I'd like to show you a more formal derivation. Note that the first formula for arbitrary (non-periodic) signals could be rewritten as

$$P_x=\lim_{M\rightarrow\infty}\frac{1}{(2M+1)N}\sum_{n=-MN}^{(M+1)N-1}|x[n]|^2\tag{1}$$

for some integer $N\ge 1$. For $N=1$, Eq. $(1)$ is identical with the first formula in your question. If you choose $N>1$ you simply add more chunks of length $N$ (instead of length $1$) to your sum if you increase $M$.

Now assume that the signal $x[n]$ is periodic with period $N$. This means that the sum over $2M+1$ periods equals $(2M+1)$ times the sum over a single period:

$$\sum_{n=-MN}^{(M+1)N-1}|x[n]|^2=(2M+1)\sum_{n=0}^{N-1}|x[n]|^2\tag{2}$$

Plugging $(2)$ into $(1)$ gives

$$P_x=\lim_{M\rightarrow\infty}\frac{1}{(2M+1)N}(2M+1)\sum_{n=0}^{N-1}|x[n]|^2=\frac{1}{N}\sum_{n=0}^{N-1}|x[n]|^2\tag{3}$$


Frequency response creates non-causal impulse response I have frequency response data. It is output from an acoustic simulator.

I convert this into an impulse response, and get:

In that picture the blue line is the impulse response. Mostly it is hiding behind the green line (which is the impulse response if I halve the frequency resolution).

I'm sorry for an unclear picture. I will replace it in due course.

Anyway, the thing I notice of interest is that if the impulse response is regarded as a loop, the ends appear to match up.

Is there any mathematical basis for this?

It seems strange to me that the impulse response is somehow non-causal.  Because what it is modelling is a single pressure wave hitting a microphone.  Commonsense suggests to me that it should have a very sharp attack as the main wave hits followed by diminishing ripples corresponding to echoes.

So I'm quite surprised to see ripples before the main spike, let alone before t=0, as it were.

Is this some artefact of the FFT due to the fact that it is forced to sample at discrete intervals, therefore the spike cannot be completely accurate, and maybe the preceding wobbles are compensating for this.

In which case would I do better to reposition the start of the impulse response 100 samples in from the end, and have it wrap round?

I am familiar with similar phenomena. During measurements with usage of sweep-sine, when you convolve your signal with inverse filter you are obtaining your IR. Although if you consider perfect case of convolving your original sweep $x(t)$ with the inverse filter $f(t)$, then you should get perfect Kronecker delta. Unfortunately what you will get is an impulse with pre and post-ringing:

One way to minimise this effect is to use so called Kirkeby Algorithm, described here by A.Farina (my guru of acoustic measurements): Advancements in Impulse Response Measurements by Sine Sweeps. I am sure will find it quite fascinating.

One last tittle-tattle:

  Commonsense suggests to me that it should have a very sharp attack as the main wave hits followed by diminishing ripples corresponding to echoes.

You are definitely right, although in real world measurements it is possible to have preceding spike before main impulse in some cases. For example if you place your source close to the wall, then wave from speaker will excite mechanical wave in a wall and wave will propagate in the material with speed of sound much higher than in case of air. Nextly wall will emit acoustic wave that will be captured by your microphone.


Converting frequency from $\textrm{Hz}$ to radians-per-sample In MATLAB I have to pass cut-off frequency for designing a filter. But this Cut-off frequency is in radians-per-sample. How do I convert my analog Cut off frequency in $\textrm{Hz}$, into the required radians-per-sample for MATLAB?

Problems like these are best attacked using some dimensional analysis:

$$f_{[\rm rad/samples]} = f_{[\rm cycles/sec]}\cdot \frac{\text{sec}}{\text{samples}}\cdot \frac{\text{rad}}{\text{cycle}}$$

$$f_{[\rm rad/samples]} = f_{[\rm cycles/sec]}\cdot \frac{2\pi}{f_s}$$

where $f_s$ is the sample rate in $\textrm{Hz}$.


multiplying the received signal by carrier, in OFDM system Sorry for this question, I'm still new in signal processing area.

In OFDM system, let $x$ be the received signal, what do we call the process of multiplying the received signal $x$ by $e^{-j2\pi F_c t}$ , where $F_c$ is frequency carrier and $t$ is interval 0 : sampling rate: sampling rate * length of signal.

So the question, what is that function? and why do we use it.

That is just a simulation of a direct conversion downmixer.

You're probably looking at a simulation of OFDM in e.g. Matlab. In a real OFDM system, the signal at carrier frequency is analog, and has to be mixed down to complex baseband with an analog mixer. Of course, in an analog mixer, there's no discrete time $t$, but just an oscillator of the right frequency with which the I and Q output branches of the mixer are generated.

Since complex baseband is so fundamentally important to how OFDM works, I strongly recommend refreshing your knowledge of it and of direct conversion mixers, if this sounds new to you.


DAC and ADC architecture in SDRs What are the DAC and ADC architectures in SDRs for inphase and quadrature phase channels?

Do we have separate ADC and DAC for each channel or there is some clever way that they are  multiplexed to reduce area.

Usually in diagrams one DAC and one ADC is shown. But then I wonder how the inphase and quadrature phase information will be separated. Would appreciate if we can get a architectural diagram with an explanation

To add here are two diagrams showing common transceiver architectures: (1) a super-heterodyne where the down-conversion is done first to an IF frequency and then to baseband and (2) a zero-IF receiver where RF is translated directly to baseband. Note in both architectures it is arbitrary (technical / technology choice) where the ADC/DAC boundary is as either approaches could be done with baseband, IF or direct RF sampling to the extent technology allows. Ultimately, at baseband we desire a complex IQ signal to support most modern modulations which require asymmetrical spectrums about the RF carrier for maximum spectral efficiency (and hence a complex baseband signal with I and Q components). We can separate into I and Q in the analog using quadrature local oscillators and then use a dual ADC/DAC, or we can use a single ADC/DAC and then separate into I and Q digitally using a quadrature NCO.

The transmit and receive architectures need not match (Can use a ZIF transmitter and Super-het receiver for example).

For high efficiency with a digital IF/RF approach, quadrature sampling can be implemented at a $f_s/4$ IF frequency with a single ADC datapath interleaving the I and Q channels as as follows:

Channel I: 1 0 -1 0 1 0 -1 0

Channel Q: 0 1 0 -1 0 1 0 -1

This concept is further illustrated in the figure below. Sampling at Fs/4:

A subsequent phase rotator would be required to remove carrier offsets, or the ADC clock itself can be in the carrier tracking loop.

Placing the digital IF at $F_s/4$ or $N F_s \pm F_s/4$ also has the advantage of simplifying ADC anti-alias filter design as the images will all be equidistant.


Representing a continuous LTI system as a discrete one I am aware that there are different ways to represent a continuous time system in discrete domain (e.g. bilinear transform, impulse invariance transform).

But my problem is as follows: Given an arbitrary LTI system $h_c(t)$, its output is given as $y(t) = (x * h_c)(t)$.

Now suppose $x(t)$ is given as:

$$
x(t) = \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc}\left(\frac{t-n T}{T}\right)
$$

i.e., Whittaker interpolation formula. I can then write a discrete-time system:

$$
y[n] = \sum_{q=-\infty}^{\infty} h[q] x[n-q] ,
$$

such that $y[n]$ is identical to $y(n T)$.

My question: What is the relation between $h[n]$ and $h_c(t)$?

There must be an exact one and it is neither given by impulse invariance nor bilinear transform.
I am aware that when $H_c(f)$ is bandlimited to $1/(2B)$, the answer is just impulse invariance transform, i.e., $h[n] = h_c(nT)$.
If this is not satisfied, it will result in aliasing but this should not prevent such a relationship from existing.

The sequence $h[n]$ consists of samples of a low pass filtered version of $h_c(t)$. Define that low pass filtered impulse response as

$$\tilde{h}_c(t)=\int_{-\infty}^{\infty}h_c(\tau)\frac{\sin(\pi(t-\tau)/T)}{\pi(t-\tau)/T}d\tau\tag{1}$$

and

$$h[n]=\tilde{h}_c(nT)\tag{2}$$

The output signal is

$$\begin{align}y(t)&=(x\star h_c)(t)\\&=\int_{-\infty}^{\infty}x(t-\tau)h_c(\tau)d\tau\\&=\int_{-\infty}^{\infty}\sum_{k=-\infty}^{\infty}x[k]\frac{\sin(\pi(t-kT-\tau)/T)}{\pi(t-kT-\tau)/T}h_c(\tau)d\tau\\&=\sum_{k=-\infty}^{\infty}x[k]\int_{-\infty}^{\infty}h_c(\tau)\frac{\sin(\pi(t-kT-\tau)/T)}{\pi(t-kT-\tau)/T}d\tau\\&=\sum_{k=-\infty}^{\infty}x[k]\tilde{h}_c(t-kT)\tag{3}\end{align}$$

Consequently, we have

$$\begin{align}y(nT)&=\sum_{k=-\infty}^{\infty}x[k]\tilde{h}_c((n-k)T)\\&=\sum_{k=-\infty}^{\infty}x[k]h[n-k]\end{align}\tag{4}$$

This results is also obvious from the fact that $(x\star h_c)(t)=(x\star\tilde{h}_c)(t)$ must hold, because $x(t)$ is band-limited. And since $\tilde{h}_c(t)$ is also band-limited, we can simply take samples from $x(t)$ and $\tilde{h}_c(t)$ to obtain the result.


What are best practices to compute an audio spectrogram? The spectrogram if generally defined with the squared magnitude of the fft. However, in lots of implementation, it seems that people just use the magnitude without square.

Moreover an audio signal is by convention scale between -1 and 1. This scaling often needs a supplementary step in implementations, in python language for example, which is not always do.

Finally, what are best practices to compute an audio spectrogram?
- square magnitude of the fft / magnitude of the fft ?
- Integer audio values / scaling (-1 to 1) audio values

EDIT

As the comments tell, these are questions without consequences if the aim is to plot an image of the spectrogram.

However, I would like to use the matrix of the spectrogram as an entry point for sound analysis and recognition. In this case, the computation process matters and I find curious that the implementations differ so often.

I think, perhaps, that the statement should be reversed: that a lot of implementations "just" use the square magnitude. This is because the DFT of a real-valued signal is typically complex, so that the square magnitude is just the product of complex conjugates, e.g.:

  X*(f) · X(f)

This requires less computation than the magnitude, which would follow the same algorithm, but with a subsequent square-root.

The square magnitude of a DFT is formally known as the Power Spectrum.


MUSIC algorithm terminology Looking for clarification on the notation described in the MUSIC algorithm in Ralph Schmidt's IEEE paper$^{[1]}$. The data model is:
$$X = AF + W$$
Schmidt defines the following:
X = received signal for each array element
A = steering vector
F = received signal
W = noise

M = number of array elements
D = number linear combination of incident wavefronts

X is dimension M x 1
A is dimension M x D
F is dimension D x 1
W is dimension M x 1

Questions

Where is the number of receivers accounted for? A practical system may have a $360^\circ$ element antenna, but only $K$ receivers, where $K \ll M$ (e.g.,$K=8$)
Where is the length of the data incorporated? Assuming $N$ samples for the received signal vector. Intuition leads me to believe that $F$ is of dimension $D\times N$
Does $F$ represent the received signal? If so, what does $X$ represent? Received signal at each element in the array?

$[1]\ $ R. Schmidt, "Multiple emitter location and signal parameter estimation," in IEEE Transactions on Antennas and Propagation, vol. 34, no. 3, pp. 276-280, March 1986

Yeah, notation is not ideal.

It is not - he assumes that each of the $M$ antenna elements is connected to its own RF chain, i.e., there are also $M$ receivers available. If you have fewer receivers you need to modify your $A$, it needs to contain the response your $K$ receivers observe given a wave from a certain direction.

He doesn't put it but yeah, $F$ should have $N$ columns (and so should $X$ and $W$) that correspond to $N$ subsequent snapshots (i.e., samples in time) of the received RF data.

The received signal, i.e., the (basedband version of the) RF data at the ouput of the $M$ antennas observed in $N$ subsequent snapshots is the $M \times N$ matrix $X$. The matrix $F$ comprises only the waveforms itself, you could think of it as the signal "before" the antenna array. But that intuition is also slightly misleading. Point is that each antenna is observing a superposition of $D$ wavefronts corresponding to $D$ far-field sources. This can be written as $A F$.

Note that this data model is a gross oversimplification of the actual physics. It ignores near-field sources/reflections, multipath, diffuse components and many things more. It still works well enough in some applications. But the physical interpretation we can put into it has its limits.


One Source Signal Separation I have a wav file which is gained by playing 2 different wav files with 2 speakers and recorded with mic. If I want to get one of the original signals using the resulting wav file and the other original wav file, is it possible to do this? I checked online and only got blind signal separation algorithm. Thank you in advance!

Let's sum your problem. Assuming the transfer function of the whole signal transmission chain from the source 1 to the recorded signal (DAC, speaker, room, microphone, ADC...) is $h_1$ (and $h_2$ from source 2 to microphone), then the signals you have complete access too are:

$(h_1 \star x_1)(n) + (h_2 \star x_2)(n)$ (the resulting signal)
$x_1(n)$ (one of the original signals).

And your problem is that you want to recover $x_2$? Then this is the same setup as this question, and certainly not a problem of source separation - since one of the signals is already known (up to a convolution). You can get an estimate of $h_1$, and recover $h_2 \star x_2$. And if you can assume that $h_1$ and $h_2$ are similar, you can try inverting $h_2$ and get a result closer to $x_2$.

It would be a true BSS problem is you had no prior information about both $x_1$ and $x_2$, in which case it would have been intractable (unless some assumptions can be made about the way these sounds are produced - eg one of them is speech and the other is background music).


Squaring signal for amplitude demodulation I read the following here. Can someone explain/elaborate on this?

  Squaring the signal effectively demodulates the input by using itself
  as its own carrier wave. This means that half the energy of the signal
  is pushed up to higher frequencies and half is shifted down toward DC

I didn't find any good source giving more information or examples of this.

If $m(t)$ is the message signal with $|m(t)|<1$, the AM signal can be written as

$$s(t)=(1+m(t))\cos(\omega_ct)\tag{1}$$

where $\omega_c$ is the carrier frequency in radians per second. Squaring $s(t)$ gives

$$s^2(t)=(1+m(t))^2\cos^2(\omega_ct)=(1+m(t))^2\frac12(1+\cos(2\omega_ct))\tag{2}$$

Low pass filtering $(2)$ gives

$$x(t)=\frac12(1+m(t))^2\tag{3}$$

from which the message $m(t)$ can be retrieved:

$$m(t)=\sqrt{2x(t)}-1\tag{4}$$


What is the meaning of channel capacity? If the transmitted information signal takes values from a modulation constellation, then is there a way to know what the channel capacity is?

Is there a connection between the modulation kind and the channel capacity?
Why and when do we do $\log_2(M)$ where $M$ is the number of modulation symbols? For example, if the modulation is 16 QAM, then what does $\log_2(16)$ tell us?
In equalization, when pilot symbols are sent, we get better estimation quality. This is known as non-blind equalization. But in many research articles often it is mentioned that non-blind technique wastes a lot of bandwidth due to periodically sending the training symbols. What is the meaning of bandwidth wastage and how does it occur? Does capacity decreases with using training/ pilot symbols?

These specific information is extremely hard to extract from text books and other online tutorials. Hence I have posted these problems here, hoping to get an intuitive answer. Shall be really obliged.

1) Is there a connection between the modulation kind and the channel capacity?

The capacity of a channel indicates the upper limit of how many bits can be transmitted per second over the channel with no errors (okay, technically it is "arbitrarily low number of errors", but it's basically the same thing).  We do various things to try to get as close to that upper limit as we can, such as modulation types, error correcting codes, etc., but none of them affect the channel capacity itself.

The equation for channel capacity is:

$$
C=B*\log_2(1+\frac{S}{N})
$$
where $B$ is the channel bandwidth in Hz, and $\frac{S}{N}$ is the signal to noise ratio (which is unitless).  Given this, the channel capacity only depends on its bandwidth and the received signal-to-noise ratio.

Going to higher order modulation types (e.g. from QPSK to 16-QAM, from 16-QAM to 256-QAM, etc.) is one way to take advantage of the $\frac{S}{N}$ portion of the equation.  If you have a very high signal-to-noise ratio it is wasteful to transmit QPSK, because the constellation points are very far apart from each other.  By going to 16-QAM the points are closer together, but they are still far enough apart to avoid most errors while transmitting twice as many bits.  If you still have "excessive" signal-to-noise you could go to a higher order QAM signal yet, allowing more bits per symbol to be transmitted.

  Why and when do we do $\log_2(M)$ where $M$ is the number of modulation symbols? For example, if the modulation is 16 QAM, then what does $\log_2(16)$ tell us?

$M$ is not the number of modulation symbols, it is the number of constellation points.  What $\log_2(M)$ tells us is the number of bits that each symbol represents.  There are 16 constellation points in 16-QAM, so we could number them from 0 to 15.  How many bits does it take to represent 0 to 15?  4 bits, or $\log_2(16)$.

  3) In equalization, when pilot symbols are send, we get better estimation quality. This is known as non-blind equalization. But in many research articles often it is mentioned that non-blind technique wastes a lot of bandwidth due to periodically sending the training symbols. What is the meaning of bandwidth wastage and how does it occur?

What they mean is that some of the bits that could have been used to transmit information are used to transmit pilot sequences instead.  In LTE, for instance, in some OFDM symbols every third subcarrier is a pilot instead of data.  If each subcarrier is using QPSK modulation, then each pilot subcarrier means that there are two fewer data bits.

  Does capacity decreases with using training/ pilot symbols?

No, not necessarily.  Yes, the number of information bits is decreased by the pilot symbols as described above, but remember that the channel capacity refers to how many error-free bits can be transmitted over the channel.  The pilots are essential for error-free communication.

A similar example is error correction codes (e.g. Turbo codes).  Error correction codes use parity bits or symbols to detect and correct errors.  We could replace those parity bits with information bits, but that would actually reduce the number of error-free bits that we sent through, so using the error correction code is a net gain.


How to approximate the transfer function of compensator using root locus? I have the root locus of a system with and without compensator. The system is consist of a PD-controller and the process, the output is fed back to the controller via $G_s(s) = 1$.
The upper graph is the root locus without compensator.
The compensator is in series with the PD-controller.
$$G_1(s) = \frac{G_r(s)G_p(s)}{1+G_r(s)G_p(s)}$$
$$G_2(s) = \frac{G_c(s)G_r(s)G_p(s)}{1+G_r(s)G_p(s)G_c(s)}$$
Here $G_r(s)$ is the transfer function of the controller, $G_c(s)$ the transfer function of the compensator and $G_p(s)$ the transfer function of the process. The system is LTI.
This is how I did it:

The second graph has two new zeros and one extra pole while contains all the zeros and poles of the first graph( system with no compensator), so the compensator looks like:
$$G_c(s) =K \frac{(s-z_1)(s-z_2)}{s-p_1}$$
But what about the constant $K$ ? I don't know how to specify it.

K is the gain parameter for which you will be at any particular point on the root locus. As you increase K, the closed loop poles will start from the open loop poles and move toward the open-loop zeros as depicted in the root locus showing all possible locations of closed loop poles versus gain K.
This is because the general form for the closed loop and open loop gain in a general loop equation is related as:
$$G_{CL}(s) = \frac{G_F(s)}{1+K G_{OL}(s)}$$

Where $G_{CL}(s)$ is the closed loop gain, $G_{F}(s)$ is the forward gain from input to output and $G_{OL}(s)$ is the open loop gain with the real gain $K$ pulled out as a independent variable that we will sweep when creating the root locus plot.
The denominator of the closed loop gain is called the "characteristic equation" and given properness of realizable causal systems (properness means the number of zeros is always less than or equal to the number of poles, not including poles and zeros at infinity), we can determine all loop characteristics (stability, dynamics, etc) from the characteristic equation alone. To really see what is occurring in the root locus plot, consider this equation expressed in terms of poles and zeros of the open loop gain $G_{OL}(s)$ as follows:
$$1+KG_{OL}(s) = 1 + K \frac{G_z(s)}{G_p(s)}$$
Where $G_z(s) = (s-z_1)(s-z_2)...$ is the polynomial of open loop zeros, and $G_p(s) = (s-p_1)(s-p_2)...$ is the polynomial of open loop poles. This is rearranged to be:
$$= \frac{G_p(s)+KG_z(s)}{G_p(s)}$$
Since this equation is in the denominator of the closed loop gain $G_{CL}(s)$ then the roots of the numerator (the zeros of this equation) will be the poles in the closed loop gain!
Knowing that, consider when $K$ is varied from very large to very small: When $K$ is very small then $G_p(s)$ dominates the equation $G_p(s)+KG_z(s)$ and therefore dominates the poles of the closed loop gain. Specifically in this case of very small $K$, the poles of the open loop gain will be the poles of the closed loop gain.
Similarly, when $K$ is very large then $G_z(s)$ dominates the equation $G_p(s)+KG_z(s)$ and therefore dominates the poles of the closed loop gain. Specifically in this case of very large $K$, the zeros of the open loop gain will be the poles of the closed loop gain.
This is why the poles start at the open loop poles and move toward the open-loop zeros in the root locus.


Digital filter not removing noise at specific frequency Matlab Doing some work at the minute on digital filters in matlab, I have a file with artifical noise added (sine wave added at specific frequency). The goal is to filter the signal and get it as close as possible to the clean signal provided.

I've done an FFT and plotted the results and found a very large spike at 29.3Hz which is not present in the clean signal.

I've tried using a notch filter, which I thought would work since it operates at such a specific frequency, however it just seems to attenuate the signal and remove some power but not block it completely. I then added a bandstop filter to try and block any signals in that region and it simply attenuated the signal also. Does anyone have any thoughts? I just seem to be lowering power of the entire signal and not actually removing anything, getting the basic shape of the clean signal but still a lot of noise present after both filters. Thanks!

[b1,a1] = iirnotch((29.3*(2/fs)),0.99999);
IIR1 = filter(b1,a1,ecg58_DC_removed);

FFT_resultFilter1 = (1/length(t))*fft(IIR1);
f=(0:1024)/1024*(200/2);
figure(4)
stem(f, 2*abs(FFT_resultFilter1(1:1025)));
xlabel ('Frequency (Hz)');
ylabel ('Spectral Magnitude');
title('First filter')
grid on

[a2,b2] = butter(2,[29.2 29.4]*2/fs, 'stop');
IIR2 = filter(a2,b2,IIR1);

FFT_resultFilter2 = (1/length(t))*fft(IIR2);
f=(0:1024)/1024*(200/2);
figure(5)
stem(f, 2*abs(FFT_resultFilter2(1:1025)));
xlabel ('Frequency (Hz)');
ylabel ('Spectral Magnitude');
title('First filter')
grid on

figure (6)
plot(t(1:1000), IIR2(1:1000));
xlabel('time (s)')
ylabel('amplitude')
title('two filters');

b1 =

1.0e-04 *

0.1571   -0.1902    0.1571

a1 =

1.0000   -0.0000   -1.0000

As mentioned in my comment, the filter returned by iirnotch is useless. From your filter coefficients you can see that the filter is only marginally stable due to two poles on the unit circle at DC and at Nyquist. Furthermore, even though the filter has a notch, it also attenuates all other frequencies quite strongly (apart from DC and Nyquist). The reason for that behavior is the extremely large bandwidth in your specification.

The figure below shows the magnitude responses of the filter you designed (top) and of a notch filter with a bandwidth BW = w0/35 (bottom) (note that the extremely large values very close to DC and Nyquist due to the poles are not shown in the top figure):

In any case, the bottom figure is what a notch filter should look like. If you tried that filter and it didn't do what you expected it to do, then the reason might be that your estimation of the noise frequency is wrong. Could it be that you got it wrong by a factor of $2$ (i.e., it would be a 60Hz hum)? [Also, doesn't the file name ecg58... suggest a disturbance at $2\cdot 29=58$Hz?]

So there might be several problems in your approach, but one is definitely the design of the notch filter, and if I may guess I would say that the other is the estimation of the noise frequency.


Discrete vs Continuous Signal Why do the same signals in the continuous-time and discrete-time differ in their frequency?

For example, if $x(t)$ is a continuous-time signal taking values only at $t= 1,2,3,4$ and $x(n)$ takes the same values as $x(t)$ at $n =1,2,3,4$.

A continuous time signal can be evaluated arbitrarily at any time,at least in theory, and can take on values at arbitrary points between any evaluation point. The Fourier transform, while not always well defined, can take on values at frequencies over $-\infty$ to $\infty$ and that value can be zero, and the signal in continuous time also takes on values over the infinite interval and those values can also be zero.

A discrete time signal only takes on values at discrete time intervals. It is only meaningful to consider values intermediate in time, intermediate to those  samples,  if for example the discrete time signal was sampled from a continuous signal, or if you have a reason (justification). In general, the values intermediate to the samples are not defined. You can consider the price of a stock or pork bellies on days where there are no trades. No one is buying or selling. The price isn’t the last traded price and it isn’t zero. It does have a Fourier Transform, computed as an infinite summation, and the frequency is uniquely (continuously defined) limited over $-\pi$ to $\pi$.

So there is only a subset of discrete time signals that we can meaningfully compare their Fourier Transforms, and those are signals that are Nyqvist Sampled.  Band limiting and periodic sampling alter the original continuous time signal, which alters the frequency spectrum, which almost answers your question.

The question of whether something is meaningful at intermediate points between samples in both time and frequency is more subtle.

I can have a periodic discrete time sequence and the frequencies in the Fourier transform are likewise values that are discrete in frequency, but are the intermediate values in frequency zero?

The most used DSP text book is  Oppenheim and Schaefer and the most recent edition has 3 interpretations of the DFT which I interpret as yes , no, and maybe.

I can also interpolate pork belly prices, but is it meaningful?

I think the answer lays with how you want to represent something.

If you want to describe a discrete time signal in the continuous time domain, a sequence of Dirac Delta functions is suitable where intermediate values are zero is OK. In other contexts, it may be meaningless.

Another example worth considering is a realization of a Poisson process.  You have a set of "arrival times" that you can store on a computer, limited by the amount of memory and numerical precision available on the machine. It is a list of numbers and the intermediate values denote no "arrivals" at those intermediate times.  The Poisson process is continuous time.  One typically doesn't perform Fourier analysis on a Poisson Process. There is some discretization in the computer representation, so one can quibble that it isn't strictly continuous but it is a practical representation of a continuous signal.  This is also an example of a continuous signal where Nyquist sampling is irrelevant.


Finding inverse z transform for two sided ROC? I have a z transform   $Y(z)=(z^2-z) /(z^2+1.3z+0.3)$

In this case i have two poles ,one at -0.3 and one at -1.0

I want to find inverse z transform of Y(z)

Following is my matlab code:

clc
clear all
close all
syms z n
num=[1 -1 0]
den=[1 1.3 0.3]
[r,p,k]=residuez(num,den)

We know that ztransform of δ(n) is 1 and ztransform of u[n] is 1/(1-z^-1)

Keeping value of r,p and k in mind,our partial fraction becomes:
$$Y(z)= \frac{2.86}{z-(-1)} -\frac{1.86}{z-(-0.3)}$$

keeping in view that ztransform of u[n] is 1/(1-z^-1)

$y[n]=2.86(-1)^nu(n) - 1.86(-0.3)^nu(n)$
This above line is inverse ztransform for ROC |z|>1

$y[n]=-2.86(-1)^nu(-n-1)+1.86(-0.3)^nu(-n-1)$

This above line is inverse ztransform for ROC |z|<0.3

My question is that,how can i find inverse z transform for ROC  0.3<|z|<1

I have also attached a snapshot of a table of common z transform and i have higlighted two cases that i have used

Since you have two poles, you have to decompose $0.3 <|z|<1$ into two ROCs, $R_A$ and $R_B$ such that: $$R_A \cap R_B \subseteq \{0.3 <|z| < 1\}$$

Your poles are $z=-0.3$, $z=-1$, so you have two choices for each ROC:

$|z| > 1$
$|z| < 1$

and

$|z| > 0.3$
$|z| < 0.3$

The requested ROC is given by $R_A \cap R_B = \{|z| < 1\}\cap\{|z| > 0.3\} = \{0.3 < |z| < 1\}$

In this case (and using your table),

$$y[]=-2.86(−1)^{}[--1]−1.86(−0.3)^{}[]$$


How to Filter Human Speech from an Audio clip? I'm looking for an open source library to detect Human speech from an audio clip, such as a wav file.

This is not a complete answer, but detection of words from the speech itself is almost always dictionary based - hence one needs to know what words would sound like. If you don't know language, you don't quite know the relation to the sounds. Hence, speech/audio recognition don't quite work.

There can be many primary features to classify:
Reference:
Automatic Language Identification: A Review/Tutorial  by Yeshwant K. Muthusamy , Etienne Barnard , Ronald A. Cole Signal Processing Magazine, IEEE Oct. 1994 Volume: 11,   Issue: 4 pp. 33 - 41

Or related paper of same author.

Phonetics:
Here you should try to look for is that you should try to collect and identify sounds based on primary phonetics. For example - you can use something like International Phonetic Alphabet
Prosodic:
Languages like English, French, German might have very close phonetics. In such case, you can also use Prosodic information which is essentially duration of sound, emphasis or accentuations, pitch and pitch variations.
Phonotactics:
Phonotactics refers to the rules that govern the combinations of the di erent
phones in a language. There is a wide variance in phonotactic rules across languages. For example, the phone cluster /sr/ is very common in the Dravidian language Tamil, whereas it is not a legal cluster in English.
Vocabulary:
If you are successful at this, you can probably make another classification of which phonetic alphabet orders are most likely in particular language. You probably don't have to recognize each word. But based on most likely phonetic symbol sequences, you can probably do language detection.


What is the threshold to prefer polyphase channelizer over FIR filtering? A wideband signal need to be processed to extract several narrowband channels from it. There are two straightforward approaches to do so:

For each narrowband channel shift frequency and filter it with FIR filter.
Implement polyphase channelizer and select only those chunks that correspond with target narrowband channels.

Considering the following extreme situations:

Need to select 1 narrowband (25 kHz) channel from a wideband channel (56 MHz)
Need to select 2048 narrowband (25 kHz) channels from a wideband channel (56MHz)

It seems to be obvious that in the first case time-domain filtration is more efficient and polyphase channelizer is more efficient in the second case.

The question is what will be more be more efficient (require fewer FPGA resources) if only 12 (25 kHz) arbitrary placed channels need to be extracted from a 56 MHz wideband signal?

A polyphase channelizer is not a special kind of filter. It is a structure that works well when using filters in multi rate settings. Polyphase is a sampling rate conversion method that leads to efficient implementations that are useful for building filter banks. The efficiency comes from only having to design one filter. The downside is that the extracted channels have to be uniformly allocated and the signal sample rate has to be an integer multiple of the channel's bandwidth.
Use a channelizer for 2 more channels that are uniformly allocated.


How is the energy of $x_1\cdot x_2$ related to the energies of $x_1$ and $x_2$? Let's say
first signal x1 = [1 2 3 4],
second signal x2 = [0.08 0.77 0.77 0.08] (Hamming window),
third signal x3 = x1.*x2 = [0.08 1.54 2.31 0.32].

Energy of x1: E1=30,
energy of x2: E2=1.1986,
energy of x3: E3=7.8165.

What is energy relation between x1,x2 and x3?

Knowing the energies of $x_1$ and $x_2$ is not sufficient for determining the energy of $x_3=x_1x_2$. What you can do is determine an upper bound for the energy of $x_3$ given the energies of $x_1$ and $x_2$ and their maximum values:

$$E_3=\sum_{k}\big|x_1[k]x_2[k]\big|^2\le\begin{cases}\max_k\big|x_1[k]\big|^2\sum_k\big|x_2[k]\big|^2=\max_k\big|x_1[k]\big|^2E_2\\\max_k\big|x_2[k]\big|^2\sum_k\big|x_1[k]\big|^2=\max_k\big|x_2[k]\big|^2E_1
\end{cases}\tag{1}$$

From $(1)$ an upper bound for $E_3$ is

$$E_3\le\min\left\{\max_k\big|x_1[k]\big|^2E_2,\max_k\big|x_2[k]\big|^2E_1\right\}\tag{2}$$

Unfortunately, this bound is generally not very tight. For your example you get

$$E_3\le\min\{19.178,17.787\}=17.787$$

but the actual value is $E_3=7.8165$.


Is maximizing hamming distance the same as minimizing correlation? For the design of an error correcting code, I might wish to maximize the distance between the codewords $\sum a_i \oplus b_i$.

For spreading sequences I'd like to minimize the cross-correlation at $t=0$ of $\sum \sigma(a_i)\sigma(b_i)$ where $\sigma(\cdot)$ maps binary $0\to1$ and $1 \to (-1)$ (i.e. what you'd get out of a BPSK matched filter when $E_b$ is normalized and there is no noise).

Are these criteria the same? It seems to me that since $\sigma$ acts as an isomorphism that turns cyclic$\pmod 2$ addition into cyclic complex multiplication, they are equivalent statements. Can you verify this and let me know where I can read more about it?

Two codewords $c_1$ and $c_2$ of length $n$, with elements in $\lbrace +1, -1 \rbrace$, and Hamming distance $d$, have a cross-correlation given by $$(n-d) -d = n-2d.$$ The reason is that there are $n-d$ bits that are equal and their product is $1$, and $d$ bits that are different and their product is $-1$.

Note that:

The larger the distance $d$, the smaller the correlation.
The correlation by itself does not tell you anything about the distance, since you also need to know $n$. For example, $1,1,1$ and $-1,-1,-1$ have correlation $-3$ but $1,1,1,1,1,1$ and $1,1,1,-1,-1,-1$ have correlation $0$.


What is difference between repeated poles and distinct poles? What is difference between repeated poles and distinct poles?
As far as i am able to understand is that repeated poles are those that have same value of both x and y coordinates while distinct poles are those that have either same value of x coordinate or y coordinate but neither both

Distinct poles do not need to share an x or y coordinate. They are classified as distinct as long as they do not share the same x AND y space. So any two poles are either distinct or repeated.

The effect repeated poles have on the impulse response of a filter is a little complicated, but the short answer is that it does change. This brief article describes the effects of repeated poles as well as providing proofs:

https://www.dsprelated.com/freebooks/filters/Repeated_Poles.html

Hope that helped.


Fourier transform of given signal This is the signal whose FT i need to find, at first i thought that i could solve this as a convolution of two rectangular pulses, but i could not find pulses that fit into this (it turns out that parts from -T/2 to 0 and T to 3T/2 are always somehow narrow). I could also try by definition but that's just too much work and very time consuming, the third option is derivative theorem, but i don't know how to use it, since i never actually used it. Any help appreciated!

HINT:

It can be solved by convolution. Take two rectangles, one that is non-zero in $[-T/2,T]$ and the other being non-zero in $[0,T/2]$. Now you just have to choose the amplitudes right and the convolution of the two signals will look like the one in your question.


What does the uncertainty principle say about recursive filters? The uncertainty principle is usually stated as a relationship between a continuous signal and that signal's Fourier transform, and says that
$$
\int_{-\infty}^{\infty} \! x^2 f(x) \ \mathrm{dx} \int_{-\infty}^{\infty} \! \xi^2 \hat{f}(\xi) \ \mathrm{d\xi} \ge \frac{1}{16\pi^2}.
$$

Apparently it is somewhat difficult to define an equivalent relationship for discrete time signals, but, for "reasonable" bandlimited signals with "reasonable" transforms, it is basically correct to say that the above holds true of the continuous signal reconstructed from the samples.  For a good overview I found the following non-paywalled paper useful:
Venkatesh, Kumar Raja, Vidyasagar: On the uncertainty inequality as applied to discrete signals, Int'l J. Math. and Math. Sci., 2006:48185, 2006. doi:10.1155/IJMMS/2006/48185.

So I also understand what this says about FIR filters: I can treat my kernel as a signal, apply the uncertainty principle and see that if I want a very narrowband filter, I'm going to have to have a lot of taps.

My question is what does the uncertainty principle say about an IIR filter?  If I have a very narrow band IIR filter it would seem that the uncertainty principle would say that I would have to trade something off.  What?  I can't figure out what the appropriate signal is that I should take the variance of.  Is it continuous reconstruction of the impulse response of my filter?  (I.e., so I should expect a very narrow-band IIR filter to have an extremely long startup transient?)

The uncertainty principle basically says the same thing for IIR filters: "the steeper the filter, the longer the Impulse Response". The only difference between FIR and IIR is that you have to use a somewhat different definition of "filter length".

While all IIR filters are technically infinite, they have very different decay rates. You can define a quantity that basically describes "energy decay over time". For most IIR filters this is an exponential decay (often dominated by the pole-pair closest to the unit circle) and the decay rate is directly related to the filter bandwidth.

An interesting quantitative approach to determining decay rates can be found in the area of reverb time calculation, using "time reversed impulse response integration".


Anti-Aliasing when oversampled From a comment in another question I read a statement

  "In addition, the initial rate of 48kHz is already performing 6 times
  oversampling and therefore the subsequent decimation by 6 will not
  cause any aliasing either."

Could you please explain it? If I would oversample by 2 and then downsample it by 2 without AA-filter I would have stuff folding back to by final Nyquist band. Is this not exactly what is called aliasing?

So the point is, oversampling with a factor of $N$ requires your signal to be confined to at most $\frac1N$ of your Nyquist bandwidth. Otherwise, you're not oversampling.

If the signal is six time oversampled, then there's no loss of information/aliasing when reducing the sample rate by six.

  If I would oversample by 2 and then downsample it by 2 without AA-filter I would have stuff folding back to by final Nyquist band.

No, you wouldn't; there's, by definition of oversampling, no energy in the "upper half" of your spectrum.

Now, through discussion in the comments it became clear that "oversampled signal" might be ambiguous:

If you had your signal alone, AA filtering would do nothing, because there's no spectral component that could be folded back (aliased) into the decimated Nyquist band. However, if we're observing noisy signals, you'd have noise energy (and be it just the quantization noise) that is relatively wideband. Thus, in that situation, although the signal of interest is oversampled, the noise is not, and to avoid getting noise aliased into our band of interest, proper AA filtering is necessary.

This is really the difference between mathematically perfect signals and noisy signals; for example, if you generated a digital sine with 100 Samples per period, you could, without incurring any problems, decimate strongly – because a digitally generated signal typically has zero noise.


What is it called when you have a single-pole recursive digital lowpass filter where the coefficients don't add up to unity? Is there a name for a single-pole recursive digital lowpass filter where the coefficients don't quite add up to 1?

For example: with the filter y[n]=a0*x[n]+b1*y[n-1] instead of using coefficients

a0=0.1
b1=0.9

What if I used

a0=0.11
b1=0.9

Is there a name for this, or is it just functionally equivalent to something else like y[n]=a0*x[n]+b1*y[n-1]+c0

Also, at what point will the filter become unstable? The example above where the coefficients add up to 1.01 looks stable to me but I think a0=0.009 and b1=1.001 would be unstable. Where's the line where it crosses from stable to unstable? b1=1?

Stability only depends on $b_1$, and it needs to satisfy $|b_1|<1$ for stability of the recursion.


What is the time reversal convolution Given a signal $r(t)$ which is a result of convolution between signal $x(t)$ and a channel $h(t)$ as below :

$r(t) = h(t)*x(t); $

what I know, the time reversal convolution can be process as follows :
$y(t) = r(t)*h^*(-t) = h(t)*x(t)*h^H(-t);$ where * denote to the convolution and $h^H$ is the conjugate. I think that right and clear, but what does mean $h(-t)$.  is it all simply equals to $conv(r(t), -h(t))$ for example as below :

x = randn(1,5);
h = randn(1,3);
r = conv(h,x);
y = conv(r,-h);

Is y is correct in the above example? however I think it should be continuous where the above example is discrete signal.

I believe that h(-t) means a "time-reversed" version of h(t). Your command: 'y = conv(r,-h);' computes the convolution of 'r' and negative 'h', and you don't want that. I think you want:

y = conv(r,conj(fliplr(h)));


Berchin's FDLS arbitrary filter design algorithm I'm trying to use Berchin's FDLS methods (for those that don't know what it is, it's a way to design arbitrary magnitude and phase of an IIR Eq)

I tried different things but neither gave me what i asked:

I tried to design an allpass filter with arbitrary phase response (nothing fancy just a linear phase and a short deviation for 5 points) and it didn't worked as planned. Actually what i want to do is to counter the phase change of an IIR filter
I tried a linear phase response for a $-12\textrm{ dB}$ notch and the same way it didn't worked (in black the desired response)

I used a samplerate of $316$, $M=158$ points and a filter order of around $50$.

Do I need to use more points to increase the filter order or is it just impossible to design my filters using that algorithm?

FDLS requires a causal frequency response. Your prototype frequency response has zero phase everywhere, which is most definitely not causal.

An IIR filter order of 50 is humongous. When FDLS has too many poles and zeroes available, it "tries" to cancel excess poles with excess zeroes. Unfortunately, due to numerical limitations, the cancellation is often pretty poor. For your notch filter, try 2nd or 4th order numerator and denominator, after making the response causal.

Though (number of measurements) = (number of coefficients) = [(numerator order) + (denominator order) + 1] is mathematically sufficient, I find that FDLS likes many more measurements than that. I recommend an absolute minimum of twice as many measurements as coefficients, and prefer 5x to 10x.

Greg Berchin


How to understand FFT results of scipy.fftpack? I calculated FFT for a speech wav-file using scipy.fftpack. How do I read (understand) the return of FFT? I have read that it supposed to be like so: y[0] is 0Hz loudness, y[1] is 1Hz loundess, ... y[n] is nHz loudness ... But seems like it is not like that exactly.

Q1: What will I get when I do abs(y)? I know that we get list of complex numbers from FFT and need to square() or abs() them. But what we will have after that? Is this Decibels?

Q2: Why do we need normalize wav-data before doing FFT? What does depend on this? Before and after normalization I get different results from FFT. If I do normalization, then results of FFT are measured by hundreds, if I don't the results are measured by 1.x small values... Is this Decibels also?

# Read wav-data
fs, data = wavfile.read('eric.wav')

# this is a two channel soundtrack, I get the first track
wavdata = data.T[0]

# this is 16-bit track, b is now normalized on [-1,1)
wavdata = wavdata / (2.0**15)

Q3: What is length of the returned list from FFT? Seems like the length of the result depends on length of given sound file... But in Q1 I supposed to get list of frequencies and their loudness independently from a given source of data. For now, if I cut in half wavdata I will get twice shorter resulting list from FFT...

Complete simple code:

import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.io import wavfile

# load the data
fs, data = wavfile.read('eric.wav')

# this is a two channel soundtrack, I get the first track
a = data.T[0]

# calculate fourier transform
y = fft(a)

# show
plt.plot(abs(y), 'g')
plt.show()

Q4: What do I process results form FFT to get it in form Db vs Hz?

Wav-file could be found here: https://aacapps.com/lamp/voices Thanks.

One rarely takes an FFT of a whole audio file, or a half. Usually only a sequence of shorter segments some number of milliseconds in duration.  The length of the "list of frequencies" will be half of the FFT length (for strictly real data) and the frequencies related to each y[i] will then be proportional to the sample rate and inversely proportional to your chosen FFT length.  Thus y[1] will only represent the output of a leaky narrow-band filter centered at 1 Hz if the FFT length is equal to the sample rate.

Decibels are more related to the log of the abs() or log magnitude of the FFT results. But without a calibrated microphone and DAC, the resulting decibels will be relative to some arbitrary level or offset. Normalization merely changes the arbitrary offset.


Audio processing -- difference between block and sample-by-sample? If a system operates on a signal with a time-domain IIR or FIR, why would an acquisition system chunk the audio into powers of two? I can understand filling a buffer by a power of two for an FFT operation. Is there a difference?

I am still confused by the difference between block processing and sample by sample processing. Some clarification would be greatly appreciated.

Note that FIR filters are sometimes implemented using an FFT (overlap-add, overlap-save). In that case it makes sense to have power of 2 buffers lengths (depending on the FFT implementation). This is of course an example of block processing, where you have to wait for a whole block before the computation of the output signal can begin. The consequence is that you always have some latency. Sample by sample processing is possible with a time-domain implementation, where you get one output sample for each input sample. Note that a time-domain implementation is not necessarily sample-by-sample but could also be using block processing.


What is usual independence assumptions on adaptive filters In this paper stands:

  The derivation and analysis of NLMS rest upon the usual independence assumptions.

It has a footnote:

  The independence assumptions used in the analysis of adaptive filters
  are:

  sequences $x(k)$ and $w(k)$ are zero mean, stationary, jointly normal,
  and with finite moments
  the successive increments of tap weights are independent of one
  another; and
  the error and $x(k)$ sequences are statistically independent of one
  another

I have several problems with that:

Who and when choose this assumptions? Is it something general or it just appears in the first publication on this topic?
I think in practice it is not possible to work just with stationary inputs $x(k)$. And also input and the adaptive weights also will never be zero mean (especially during online usage when it is impossible to normalize the input data). Also the error and input are strongly correlated, so how it can be statistically independent?
Is there some other theory dealing with analyzing those more real-life applications where those assumptions cannot be fulfiled?

If is it possible to answer, than please answer. If my questions are too stupid, please correct me. Thanks in advance.

I think there is an error in your referenced independence assumption.
$w(k)$ should be the update part $\Delta w(k)$ i.e the

$w(k+1)=w(k)+\mu \Delta w(k) =w(k)+\mu \frac{x(k)*e(k)}{c+x(k)^Hx(k)}$

The error after convergence is uncorelated and zero mean.

The two other assumptions are correct.

If your filter is fast enough you may use it at nonstationary systems as long as you are able to track the variations.


Meaning of Real and Imaginary part of Fourier Transform of a signal Say $f$ is a signal of time $t$, $F$ its Fourier transform of the variable $v$.

It is known that in polar coordinate, $|F(v)|$ tells us how much the frequency $v$ is present over the signal, and $Arg(F(v))$ tells us how much the contribution of this frequency is phase-shifted.

What information do its real and imaginary part tell us ?

Or if I reformulate my question : can we give an interpretation of the Fourier transform in Cartesian coordinate like we can do in polar coordinate ?

The real and imaginary parts of the Fourier transform of a signal $x(t)$ are the Fourier transforms of the signal's even and odd parts, respectively:

$$X_R(\omega)=\frac12[X(\omega)+X^*(\omega)]\Longleftrightarrow\frac12[x(t)+x^*(-t)]=x_e(t)\\
X_I(\omega)=\frac{1}{2j}[X(\omega)-X^*(\omega)]\Longleftrightarrow\frac{1}{2j}[x(t)-x^*(-t)]=-j\cdot x_o(t)
$$

where $X_R(\omega)$ and $X_I(\omega)$ are the real and imaginary parts of $X(\omega)$, and $x_e(t)$ and $x_o(t)$ are the even and odd parts of $x(t)$, respectively.


Is filtering necessary if I'm not downsampling? My data is sampled at $f_s$, and I'm interested in analyzing certain frequencies $0<f_1<\ldots<f_N<\frac{f_s}{4}$. I then take $T$ samples of data, $\{x_1,\ldots\,x_T\}$ where $T\gg \frac{1}{f_s}$, FFT it after padding to the next power of 2 after $T$ and then pick the FFT bins closest to my $f_i$ of interest to proceed with the next step.

My question is: Is it necessary to (or is there any advantage if I) bandpass filter (BPF) my data between $f_1$ and $f_N$? I figured that since there was no downsampling involved, there won't be any aliasing and I needn't BPF it. Besides, a well designed filter should not modify the frequencies of interest, so it makes no difference. Am I missing something or is filtering a necessary step before FFT?

Note: I won't be inverting it back to the time domain. I'm just interested in the specific bins (but several that FFT is much faster than Goertzel).

Filtering is not necessary as it would only effect those frequency bins outside of the passband.

However, depending on how small your passband is, you could potentially save yourelf a lot of FFT processing by low-pass filtering and then downsampling before your perform your FFT.  For example, if you downsample by 4, the resulting FFT could be 4 times shorter and still have the same frequency resolution.


Calculating the gain for adjusting sound level to X db I have a sound signal, and I would like to implement a simple realtime level control (for info, I am using Pure Data, but I'd like to know about the general case).

For this I thought of using the envelope amplitude (in dB) of the input Vin_db, calculated over a window of N samples ; use it to calculate the gain G to apply to the signal to get an output with amplitude Vout_db. (Disclaimer : I am very rusty with maths and decibels ...)

So I took the wikipedia formula for decibels :

V_db = 10 * log10(P / P0)

=>

Vin_db = 10 * log10(Pin / P0)

Vout_db = 10 * log10(Pout / P0)

And this gave me :

Pout = Pin * 10^((Vout_db - Vin_db) / 10)

=>

G = 10^((Vout_db - Vin_db) / 10)

However, this doesn't seem to work ... When I measure the output level it is not the expected Vout. Is there a fundamental flaw in my reasoning ? Or just a calculation error ?

EDIT

Thanks to Paul R I found the correct formula (I need a 20 instead of a 10):

G = 10^((Vout_db - Vin_db) / 20)

You don't need to convert the voltages to dB, just calculate the required gain (as a factor) and then convert this to dB if necessary, i.e.

G = Vout / Vin                 # required gain (factor)

G_dB = 20 * log10(G)           # required gain (dB)

where

Vin = input voltage
Vout = output voltage

If you really want to work with dB throughout then you can do it like this (should give the same answer):

Vin_dB = 20 * log10(Vin)       # input voltage (dB)

Vout_dB = 20 * log10(Vout)     # output voltage (dB)

G_dB = Vout_dB - Vin_dB        # required gain (dB)

G = 10 ^ (G_dB / 20)           # required gain (factor)


If low frequency travels longer distance, why our speech is not travelling longer distance? It is know that low frequency signals travel longer distance that of high frequency signal. If so, our speech signal is low frequency signal, why it is not travelling further?

Well, first of all the Sound Level Pressure decreases by $6 \; \mathtt{dB}$ when doubling the distance - this plays a big role. We do also have sound attenuation coming from our medium - air. Let's take a closer look onto sound absorption coefficient for different frequencies:

Knowing that human speech is mostly concentrated at the range of $300\;\mathtt{Hz}$ - $3.3\;\mathtt{kHz}$ (yes I know - let's forget about sibilants) you can judge that we are really doing our best. Larynx tones and first formants of voiced speech are clearly transmitted with minimum lost of energy. Therefore main problem is the decrease of SPL with increasing of distance from source.

Another issue is the hearing curve of human ear which makes the low attenuation at the low frequencies worse. For example even though air absorption for $2 \;\mathtt{kHz}$ - $3 \;\mathtt{kHz}$ is rather high, we can still hear the ambulance siren as it is in frequency range where human ear has it's resonance. Obviously that depends on many factors, i.e. in city you have lots of barriers for sound.

For example let's say that you are a male speaking with normal level. Exemplary average $dBA$ values in frequency ranges are following:

$250\; \mathtt{Hz}$: $51 \; \mathtt{dBA}$
$1\;\mathtt{kHz}$: $49\;\mathtt{dBA}$

Assume that these are values measured at the distance of $1\mathtt{m}$ We want to know how loud you are at $100\mathtt{m}$. Applying atmospheric attenuation factor for 10% of relative humidity gives:

$L_{250}^{a}=51-0.1=50.9\;[\mathtt{dBA}]$

$L_{1000}^{a}=49-1.4=47.6\;[\mathtt{dBA}]$

Not a big change, isn't it? Now let's apply the following formula for calculating the sound level at a given distance $d_{ref}$ (assuming cylindrical wave):

$$L_d=L_{ref}-20\log{\left( \dfrac{d_{ref}}{d} \right) } $$

Applying it to previously calculated values and substituting $d=1\;\mathtt{m}$, $d_{ref}=100\;\mathtt{m}$ we get:

$L_{250}^{a,d}=L_{250}^{a}-40=10.9\;[\mathtt{dBA}]$

$L_{1000}^{a,d}=L_{1000}^{a}-40=7.6\;[\mathtt{dBA}]$

Both values we've obtained are way below the level of human breathing at a distance of $1\;\mathtt{m}$  ($25\; \mathtt{dBA}$) and environmental noise. You can also try to shout but it still won't give you  big increase in a distance. Please refer to this article: Average Speech Levels and Spectra in Various Speaking/Listening Conditions for detailed info about levels of speech.

One last word, bit out of scope of this question. There are some cases when you can be heard in far distance or you could be heard very poorly in a very close distance. Such phenomena occurs a lot in water (but also in the air), where sound waves can travel tens of kilometres. This is due to presence of wind, temperature inversions, terrain profile and refraction of sound. If you wish to know more, then I am sure there are plenty of guys in Physics SE willing to answer that question. As for start, here is some interesting article: Outdoor Sound Propagation in the U.S. Civil War, as well as bunch of papers regarding Outdoor Sound Propagation.


Advancing a causal IIR filter's impulse response Given an infinite impulse response (IIR) filter as coefficients in a topology of your choice, what is a straightforward way to modify the filter so that the impulse response is shifted in time (advanced) by one sample to the left, discarding the sample that would make the filter non-causal?

Or is it necessary to modify the impulse response and solve the new coefficients from scratch?

(I asked this question because I wanted to know if it's possible to snip away, with a time shift, an arbitrary number of the initial samples of an IIR filter's impulse response without increasing the filter's complexity.)

Let's consider a causal IIR filter with the following transfer function:

$$H(z)=\frac{B(z)}{A(z)}=\frac{\displaystyle\sum_{n=0}^{M}b[n]z^{-n}}{1+\displaystyle\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

Note that $M$ need not be equal to $N$. Let $K$ be the number of samples that must be removed from the impulse response corresponding to the system $(1)$, so the new desired impulse response is

$$\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots\tag{2}$$

Let $H_K(z)$ denote the FIR transfer function corresponding to the first $K$ samples of $h[n]$:

$$H_K(z)=h[0]+h[1]z^{-1}+\ldots +h[K-1]z^{-(K-1)}\tag{3}$$

The transfer function of the new filter is then given by

$$\tilde{H}(z)=z^K\left[\frac{B(z)}{A(z)}-H_K(z)\right]=\frac{z^K\big[B(z)-H_K(z)A(z)\big]}{A(z)}\tag{4}$$

As expected, the denominator polynomial remains unchanged, we just get a new numerator polynomial

$$\tilde{B}(z)=z^K\big[B(z)-H_K(z)A(z)\big]\tag{5}$$

Let's first consider the case $K=1$. In this case, the FIR transfer function $H_1(z)$ defined in $(3)$ is trivially given by

$$H_1(z)=h[0]=b[0]\tag{6}$$

Consequently, the new numerator polynomial becomes

$$\tilde{B}(z)=\big(b[1]-b[0]a[1]\big)+\big(b[2]-b[0]a[2]\big)z^{-1}+\ldots +\\+\big(b[N]-b[0]a[N]\big)z^{-(N-1)}\tag{7}$$

if $M=N$. If $M>N$, there are additional terms $b[k]z^{-(k-1)}$, $k=N+1,\ldots,M$, and if $M<N$, there are additional terms $-b[0]a[k]z^{-(k-1)}$, $k=M+1,\ldots,N$.

The computation of the new numerator coefficients according to $(7)$ is very simple and efficient. For $K>1$ this procedure can be applied iteratively. In Matlab/Octave this is (almost) a one-liner:

for k = 1:K
    b = [b(2:M+1);zeros(N-M,1)] - b(1)*[a(2:N+1);zeros(M-N,1)];
    M = length(b) - 1;
end

Note that for any $K$, the order of the new numerator polynomial $\tilde{B}(z)$ equals $\max(M-K,N-1)$.

Original Answer (correct, but the presentation and solution given above is probably better; the examples below are still valid and hopefully useful)

A solution can be obtained by recomputing the numerator coefficients as follows:

$$\tilde{b}[n]=(\tilde{h}\star a)[n],\qquad n=0,1,\ldots,L-1\tag{3}$$

where $\star$ denotes convolution, $L=\max(M,N)$, and $\tilde{h}[n]$ is related to the impulse response $h[n]$ of the original IIR filter by

$$\tilde{h}[n]=h[n+1],\qquad n=0,1,\ldots,L-1\tag{4}$$

The method can be summarized as follows:

Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=1,2,\ldots,L$
Compute the new numerator coefficients from the impulse response and from the denominator coefficients using Eq. $(4)$

I chose a random filter to show the result of the suggested procedure in the figure below (blue: original impulse repsonse: red: new impulse response shifted to the right by one sample for comparison):

The left column below shows the original numerator coefficients alongside the modified numerator coefficients, aligned in such a way that it becomes obvious that for higher indices the new coefficients equal the original coefficients. (This is the case because in this example $M>N$; c.f. Eq. $(7)$ of the new answer above, including the text below Eq. $(7)$).

  -0.67010
   0.65673   0.92477
  -0.61219  -0.74621
  -0.60645  -0.53944
  -0.50316  -0.50316
   0.82855   0.82855
   0.20396   0.20396
   1.15680   1.15680
  -0.29830  -0.29830
   0.10605   0.10605

The denominator coefficients of both filters are

a = [1, 0.4, -0.2, 0.1]

EDIT:

This method can be generalized to remove $K$ initial values of the original impulse response, even if $K>M$:

Compute $L=\max(M,N)$ coefficients of the original impulse response $h[n]$ for indices $n=K,K+1,\ldots,K+L-1$
Compute the new numerator coefficients according to Eq. $(4)$ with $\tilde{h}[n]=h[n+K],\qquad n=0,1,\ldots,L-1$

The resulting numerator coefficients will generally have several trailing zeros (or values sufficiently close to zero) which can be removed.

The following example shows the result of removing $20$ initial values of the impulse response of the filter with coefficients

[b,a] =

   1.000000   1.000000
   2.000000  -0.187178
   3.000000  -0.075960
  -1.000000   0.812250

The numerator coefficients of the new filter are

   0.95236
   0.90498
   0.42252

Of course, these coefficients can no longer be interpreted as (an approximation of) a shifted version of the original numerator coefficients.


Confused about filter's cutoff frequency in rational sampling rate converter A interpolator which increase sampling rate by L times is followed by an anti-imaging filter, the filter is said to have a cutoff frequency of Fs/2L. Is Fs here is the frequency before or after interpolation?

Also during rational sampling rate conversion, the filter in between the interpolator and decimator is said to have cutoff at min(Fs/2L, Fs/2M). But clearly Fs isn't same for both of them. The interpolator changes the sampling rate before the signal passes through the decimator.

I assume you'r talking about digital interpolators.

An interpolator is basically an L times expander followed by an anti-imaging lowpass filter with gain L. The cutoff frquency of the lowpass filter is defined in the discrete-time frequency domain as

$$ w_c = \frac{\pi}{L}. $$

The analog equivalent of this cutoff frequency is

$$ f_c = F'_s \frac{ w_c }{ 2\pi} = L \times F_s \frac{ \pi / L}{ 2\pi} = F_s / 2. $$

This frequency refers to the increased sampling rate $F'_s = L \times F_s$.

A similar argument holds of the rational sample rate conversion system where an expander (by L) is followed by a lowpass filter and then followed by a decimator (by M).

The cutoff frequency of the intermediate lowpass filter is $$\omega_c = \min(\pi/L, \pi/M).$$ Then the analog equivalent of this lowpass filter will be
$$ f_c = L \times F_s \frac{ \min(\pi/L, \pi/M)}{ 2\pi} .$$


Calculating Beam pattern, Directivity of microphone array As answered in Relationship of Fourier Transform and Spatial Directivity Pattern the relation between directivity pattern and weights is given.  But I am bit confused, I have Beamformer weights (complex valued, for each frequency bin) for a particular steering vector. The relationship mentioned in the above link needs the values to be calculated for each theta, but the beamformer weights are for a particular steering angle theta.
Do I need to calculate beamformer weights for each of the theta, or I have a wrong understanding of the directivity pattern.
Suggest methods for calculating directivity using beamformer weights.
Thanks

You choose the weights the steer the main beam of the antenna in the particular direction you want - let's call this $\theta_s$ i.e. steering angle. Now, given that you are steering the antenna in a given direction you are also interested in the array's response to incoming signals from other directions. This is the directivity pattern and is a function of the angle $\theta$ of the source signal. For conventional beamforming (non-adaptive) the magnitude of the directivity pattern will reach a maximum when $\theta_s=\theta$.


Software defined radio panadapter I am trying to make a panadapter for a software defined radio and I am a little bit stuck and hope someone can help me.

I have an I Q signal from my radio. (Actually I am using an IQ wav file from the internet recorded at 44800 2 channels 16bit).

If I use HDSDR the spectrum looks like this:

The display from my app looks like this:

I am doing something wrong and I do not know exactly what (I am new to this and I am not very sure I completly understand the concept behind...) I will try to describe what I am doing and hope you can help me fix it

I open the file and start reading from it in a byte buffer (1024 bytes at a time)
I convert every 2 bytes into shorts (ByteOrder.LITTLE_ENDIAN)
I apply HanningWindow to the buffer
I split the buffer in 2 buffers each for a channel (the %2==0 as
left and %2==1 right )
I create and array of Complex numbers considering the left[i] as the
real part & right[i] as imag part
I do FFT on the complex array (wavenumber table with size 512)
I do an FFT shift (from [1,2,..,n/2,n/2+1,..,n-1,n] to [n,n-1,..,n/2+1,1,2,..,n/2] )
I plot the result

Do I need to process the IQ signal before somehow (demodulate it ?)

I am reading in a wrong way the file? (I know that the file is starting with header & everything but after that the data part starts)

How can I do a pitch /amplitude correction to remove the unwanted image existing in a not ideal IQ recorded by soundcard?

Thank you,
Bogdan

There are a number of things that could be going wrong; you haven't provided enough information to definitively diagnose your the problem. (A link to the wav file you obtained would be useful, so we can try it out ourselves.) Here are some observations:

You're assuming the file is uncompressed 16-bit little-endian samples, and ignoring the presence of the header. This could be wrong, either because the format is not what you think, or the header is not a multiple of 4 bytes long. Instead, use an audio-editing program such as Audacity to read the file and write it out again as a raw file where you have specified the parameters. That, or use a WAV-reading library in your program.
Your lengths don't match up. If you read 1024 bytes, then you have 512 shorts, or 256 complex-shorts — but you say your FFT size is 512, not 256. Either your program is not as you say or it might be reading uninitialized memory.
Your output appears to have four almost-copies of the spectrum. Two mirrored copies are explainable if you are misreading the input such that the I and Q components are not aligned properly — a symmetric spectrum is to be expected if one of I or Q is zero, constant, or otherwise not properly aligned with the input.

The other copy might be due to an error in your FFT shift code or graphics code, or because you are trying to display both components of the complex FFT output separately rather than taking the magnitude.
I would suggest adjusting the scaling of your color shading so that it is more similar to the HDSDR output, and increasing the FFT size (bin count). (You may also need to take the log of the magnitude.) This will make the result look more like HDSDR's, which will help you match up parts of the spectrum.

  Do I need to process the IQ signal before somehow (demodulate it ?)

No. The sorts of processing that could be done to refine a panadapter display would be done after the FFT. The only things to do before would be things like IQ imbalance correction, which are just as appropriate for a signal being demodulated and played as audio (or whatever sort of information) as one going to the panadapter.


Optimal time for sampling symbol in digital receivers I'm studying various publications about symbol timing recovery algorithms, Early-late, Mueller and Muller & Gardner. All documents I've read tells that, to minimize ISI :

  The optimal time for sampling the signal is at the peak of the receiver matched filter

Looking at a real signal, it seems that a local peak is not always the right sampling time. For example, this figure shows the output of the receiver matched filter and the ideal sampling times :

In this example, the ideal sampling times are not always at local peaks of the signal. Are there special requirements so that the "optimal sampling time is at local peaks" or could you point me in the right direction to understand this particular point ?

You are mixing up two different notions of sampling.

In a digital
communications system, the received (analog, continuous-time) signal
is passed through an A/D converter so that the needed further
processing (e.g. matched filtering) can be implemented on a
programmable DSP processor or as a MATLAB or C++ program or
on a special purpose ASIC. The sampling that you show in your question
is the one being carried out in the A/D converter, and as long as
the sampling is above the Nyqvist rate, the fact that the samples
are not at the peaks of the signal is irrelevant: the original
signal is represented perfectly adequately by the samples regardless
of where the peaks in the signal are with respect to the sampling
instants.

The other notion of sampling has to do with the matched filter
output which is a sequence of digital data. Which of these
data should the decision-making device use to decide whether
(in the simplest case) whether a 0 or a 1 was transmitted?
The answer is what you quote: it is the datum when the
signal output of the matched filter has a peak. That is,
we are picking one value from the matched filter output
sequence and making the decision based on that.
Note that

We are not going to average the matched filter output (or
take a weighted average of the output) and make a decision on that
quantity.
We are not going to take the maximum value of the actual output
sequence (which may, because of noise, be at a different location
than the time when the signal output of the matched filter is
supposed to peak) and make a decision based on that.

For more information on matched filters (including explanations
of the above bulleted points), see
this other answer of mine.


Convert one sound sample in different sample where each sample's pitch is a note (do-re-mi...) For example, I would like to do like in this video:
https://www.youtube.com/watch?v=YBhneZuAb5k

We have for example a goat yelling sample and from this sample, we can change the pitch according to the pressed piano key. I would like to make a script that will, from one input sample, generate all samples for each pitch/note (do re mi...).
What is the algorithm for doing that?

One generic term for the class of libraries and algorithms that do this, make a set of music keyboard notes from random (psuedo-periodic) pitched sounds, is "time pitch modification".  This is sometimes accomplished by frequency domain phase vocoder analysis-resynthesis or by time domain synchronized overlap-add, followed by filtered resampling.  One example library might be Apple's OS X AVAudioUnitTimePitch effect unit.  Here's a wikipedia page on the topic: https://en.wikipedia.org/wiki/Audio_time-scale/pitch_modification

First, you might have to use a pitch estimation algorithm to determine a perceived pitch for the original sound sample.  Here's a wikipedia page summarizing some of these methods: https://en.wikipedia.org/wiki/Pitch_detection_algorithm


For linear IIR digital filter, what happens for negative frequencies? By negative frequency, I refer to Fourier transform. Often, the frequency response of a digital filter is only displayed for positive frequencies. For a linear IIR digital filter, what happens for negative frequencies? Are frequency response for negative frequencies mirror images of what happens for positive frequencies?

It may be easier to start with analog signals. If the filter's impulse response $h(t)$ is real, then its frequency response $H(f)$ is conjugate symmetric: $H(-f)=H^*(f)$, where $(\cdot)^*$ indicates complex conjugate. Since $H(f)=|H(f)|e^{j\angle H(f)}$, then conjugate symmetry implies that $|H(f)|=|H(-f)|$ and $\angle H(f) = - \angle H(-f)$. This is true for any kind of filter, FIR or IIR.

So, if the filter input $x(t)$ is $$x(t)=A_0\cos(2\pi f_0t+\phi_0)=\frac{A_0}{2}\left(e^{j2\pi f_0t}e^{j\phi_0}+e^{-j2\pi f_0t}e^{-j\phi_0}\right),$$ then the filter output is

$\begin{align*}
y(t) &= \frac{|H(f_0)|A_0}{2}\left(e^{j2\pi f_0t}e^{j(\phi_0+\angle H(f_0))}+e^{-j2\pi f_0t}e^{-j(\phi_0+\angle H(f_0)}\right) \\
&= |H(f_0)|A_0\cos\left(2\pi f_ot+\phi_0+\angle H(f_0)\right).
\end{align*}$

As you can see, the amplitude of negative frequencies are affected exactly the same as positive frequencies. The phases are affected in a similar fashion, but with phase of opposite sign for negative frequencies. Note that this is necessary to preserve the symmetries that real signals have in the frequency domain.

In the case of discrete signals and filters, the conclusion regarding negative frequencies is the same. Finding the actual gain and phase at a particular frequency is slightly more difficult, though. The reason is that the DFT of the filter's impulse response defines the filter's gain and phase only at a finite number of frequencies. In general, the filter's gain and phase at an arbitrary frequency $f_0$ can be found by evaluating the filter's transfer function at $z=e^{j2\pi f_0}\,$.


Segmenting pills from the background I just recently started out with Image processing and taken up a course related to it in Grad school. But I already have a project to do without too much information about the subject, however I've made some steady progress. I am trying to segment the pills from their respective backgrounds. For images with contrasting backgrounds, I've been able to segment the pills using Otsu's method. With regards to images with similar backgrounds, Otsu's method doesn't quite work. I've been reading quite a few papers with regards to segmenting, but most of the papers that I've read use a manual thresholding value depending upon the kind of image. Is it possible to detect the correct thresholding value and automatically threshold an image and use techniques such as seed-growing or clustering to segment the image ?

The color space I've been using is Lab*, so I'd appreciate it if you could recommend the appropriate color space that I should be using too.

The images in question:

original 1

result 1

original 2

result 2

original 3

result 3

If you want to use the thresholding approach, you should use an adaptive thresholding method if there are big lighting variations like in the 3rd example image (dsp question here).

Also, you should experiment with colorspaces, it's easy: the script to decompose the image into different colorspaces should be no more than a few lines long, and a lot of image viewers have that option available. The best one should be easy to determine just visually. If you want to read up on colorspaces, there's anothrt good dsp question here.

Finally, you might want to try a diffrent approach. One idea would be to do a not-perfect segmentation, then an edge detection, and finally using something like Hough transform for circles which also works fine on (incomplete) circle arches. (this idea is ofcourse only applicable to round pills)


Zero-order-hold for two dimensional signal $x(t_1,t_2)$ From Wikipedia:

  A zero-order hold reconstructs the following continuous-time waveform
  from a sample sequence $x[n]$, assuming one sample per time interval
  $T$: $$x_{\mathrm{ZOH}}(t)\,= \sum_{n=-\infty}^{\infty} x[n]\cdot
 \mathrm{rect} \left(\frac{t-T/2 -nT}{T} \right) \ $$ where
  $\mathrm{rect}()$  is the rectangular function.

We can refer to $\mathrm{rect}()$ also as one-dimensional rectangular function. For example, the product of two one-dimensional rect functions can be viewed as a rect function in two dimension, i.e. a function
that has the value 1 on the square of side length 1 centered at the origin, and has the value 0 outside this
square. I am wondering if there exists, for example, a two dimensional version of ZOH for a signal in two variable, $x(t_1,t_2)$.

Thanks in advance.

Yes there exists a 2D analog of the 1D zero order hold for image reconstruction, whose mathematical formulation follows closely that of 1D case as:
$$I_{\mathrm{ZOH}}(x,y)\,= \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} x[n,m]\cdot \operatorname{rect} \left( {x-nT_x,y-mT_y} \right) $$

Where $x[n,m]$ is the spatial samples of 2D image to be reconstructed, $T_x$ and $T_y$ are sampling periods in $x$ and $y$ dimensions (which are generally equal in practice) and the $\operatorname{rect}()$ function is a 2D rectengular pulse of height one and durations $T_x$ and $T_y$ along the dimensions.

The zero order hold interpolation can easily be observed when you enlarge an image through an image editor by not using any smooth interpolation methods but just repeated filling of those new enlarged pixels.

In particular just use, for example, IrfanView to open any image and choose Image->Resize from menu, then choose to enlarge the image from the options in the coming window, and finally choose "Resize (faster, lo quality)" option from the "Size Methods" sections about lower right corner (for version 4.40)

If you instead choose "smooth" methods, they will fill in the pixels with calculated intermediate values, therefore implementing higher order holds than zero.


Is it necessary to apply some window method, to obtain the FFT? Java I read that the window is needed to apply filter (Low Pass, High Pass, Band Pass or band Stop),  but not for the FFT.

1.
What advantage does apply a window or disadvantage does not apply?

2.
What is the criteria to choose the size (My LengthW:Length of Han Window)?

3. Where I can find literature that specifically?

My code:

Suppose an array of samples stored in

integer SIZE = 16384; //SomeValue
short[] sSamples = new short[SIZE];
sSamples = obtainSamples();  //Some method to fill my array

Now I need to obtain my FFT.
I arbitrarily chose 4096 for my FFT, then are 8 times (or 8 fft's)

FFTLength = 4096;
SIZE = 16384;
double Overlap = 0.5;
SIZE/(FFTLength*Overlap) = 8;

DIRECT METHOD:

complex cSamples = new complex[SIZE];
cSamples = short2Complex(sSamples);

Now I need to obtain part of array of sample 8 arrays of length:SIZE/4

//getPartArray(from, until) maybe like System.ArrayCopy
complex[] cSamples00 = getPartArray(0,4095);
complex[] cSamples01 = getPartArray(4096,8191);
//similarly for all, from 00 until 07
complex[] cSamples01 = getPartArray(12287,16383);

Finally, Obatin the FFT

//fft(complex[], length)
complex[] cFFT00 = fft(cSamples00,4096);
//similarly for all, from 00 until 07
complex[] cFFT07 = fft(cSamples00,4096);

WINDOWED METHOD:

//int LengthW = ??  // (besides positive and odd)
double[] Wn = new double[LengthW];
//Suppose Han Window,can be Hamming, Kaiser, Blackman, etc
Wn = obtainHanWindow(LengthW);
double[] sWindowedSamples = applyingWindow(Wn,sSamples);

complex cSamples = new complex[SIZE];
cSamples = short2Complex(sWindowedSamples);

Now I need to obtain part of array of sample 8 arrays of length:SIZE/4

//getPartArray(from, until) maybe like System.ArrayCopy
complex[] cSamples00 = getPartArray(0,4095);
complex[] cSamples01 = getPartArray(4096,8191);
//similarly for all, from 00 until 07
complex[] cSamples01 = getPartArray(12287,16383);

Finally, Obatin the FFT

//fft(complex[], length)
complex[] cFFT00 = fft(cSamples00,4096);
//similarly for all, from 00 until 07
complex[] cFFT07 = fft(cSamples00,4096);

The only difference is calculate a window and apply it to the input.

Definition's class and methods

The complex class is defined like:

class complex {
  public double real = 0.0;
  public double imag = 0.0;
}

short2Complex:

short2Complex(short[] intSSample) {
  complex intCSample = new complex[intSSample.length];
  for (int i = 0; i < intSSample.length; i++) {
    complex tcomplex = new complex();
    tcomplex.real = intSSample[i];
    intCSample[i] = tcomplex;
  }
  return intCSample;
}

obtainHanWindow:

double[] obtainHanWindow(int Size) {
  //Size must be Odd and Positive integer
  double[] dHan = new double[Size];
  int HalfSize = (int)Math.floor((double)Size/2);
  for (int i = 0; i < HalfSize; i++) {
    int j = 2*HalfSize - i;
    double Pi2DivSize = 2.0*Math.PI/Size;
    dHan[i] = 0.5+0.5*Math.cos(Pi2DivSize*(double)(i - HalfSize));
    dHan[j] = 0.5+0.5*Math.cos(Pi2DivSize*(double)(HalfSize - i));
  }
  dHan[HalfSize] = 1.0;
  return dHan;
}

applyingWindow

applyingWindow(double[] dWn, double[] dSamples) {
  double dOutput = new double[dSamples.length]
  for (int i=0; i < dWn.length; i++) {
    dOutput[i] = 0;
    for (int j=0; j < dSamples.length; j++) {
      if(i-j>=0) {
        dOutput[i] += dSamples[j]*dWn[i-j];
      }
    }
  }
  return dOutput;
}

As you can see; When a sinusoidal signal of arbitrary frequency is input to the ‘N’ point DFT, the spectral response of the signal will not be a Delta function (impulse) because when we take into account of the missing end point to be the same as the beginning of the next period of the periodic extension of the sequence, the periodic extension of the signal is no more a sinusoid. One way to solve this problem is to ensure Coherent Sampling as shown in figure 1

But coherent sampling can be done only if the DFT engine receives pure monotone signals. So in a real world application, Coherent sampling is not feasible.

SO, just diminish the magnitude of the discontinuities by multiplying the whole signal with some function.

So, when you multiply by window function, it eliminates the jump discontinuities and removes the false frequency that is introduced when coherent sampling is not done (see Figure 3.

The length of the window should be the length of your signal.
Selection of your window could be depending on how much side-lobe height do you tolerate.

Image Source: http://electronicdesign.com/analog/choose-right-fft-window-function-when-evaluating-precision-adcs


How to prove $\int_{-\infty}^{\infty}x^2(t)dt=\int_{-\infty}^{\infty}\hat{x}^2(t)dt$ where $\hat{x}(t)$ is the Hilbert Transformation of $x(t)$ $\hat{x}(t)$ is the Hilbert Transform of $x(t)$. Can anybody help me proving the above.

The Hilbert transform of $x(t)$ denoted by $\hat{x}(t)$ has the following Fourier transform
\begin{equation}
 \mathcal{F}(\hat{x}) = H(w)X(w) = -j \operatorname{sgn}(w)X(w)
\end{equation}
because
\begin{equation}
 \hat{x}(t) = \int\limits_{-\infty}^{\infty} \frac{x(u)}{\pi(t-u)} \ du
\end{equation}
Notice that
\begin{equation}
 \vert \mathcal{F}(\hat{x}) \vert^2 =
 \vert -j \operatorname{sgn}(w)X(w) \vert^2
 =
 \vert X(w) \vert^2
\end{equation}
i.e.
\begin{equation}
\int\limits_{-\infty}^{\infty} \vert \mathcal{F}(\hat{x}) \vert^2
\ dw
=
\int\limits_{-\infty}^{\infty}
\vert X(w) \vert^2
\ dw
\end{equation}
Use Parsevals equation which gives you
\begin{equation}
\int\limits_{-\infty}^{\infty} \vert \hat{x}(t) \vert^2
\ dt
=
\int\limits_{-\infty}^{\infty}
\vert x(t) \vert^2
\ dt
\end{equation}
P.S: For real signals, you could indeed remove the magnitudes (absolute values).


Zero-Phase Algorithm in Matlab From this document it's says

1.Filter the Data
2.Flip that
3.Filter again
4.and Flip that again

for zero-phase filter. but when i do this in matlab , i'm compare between filtfilt method and what i did , it isnt equal. what is the problem ?

MATLAB CODE

clc;clear all;
emg_healthy=load('emg_healthy.txt');
x= 10;
fs= 4000;
my_counter = 1000;
myRawData = emg_healthy(:,2);
myRawData = myRawData(1:my_counter);
myDataRecti = abs(myRawData);
[b,a]= butter(6,x/(fs/2),'low');
filter_y=filtfilt(b,a,myDataRecti);
firstFilter = filter(b,a,myDataRecti);
flip_firstFilter = fliplr(firstFilter);
secondFilter = filter(b,a,flip_firstFilter);
flip_secondFilter = fliplr(secondFilter);
xlabel('Sample Number');ylabel('Genliği');plot(flip_secondFilter,'b');
hold on ;
xlabel('Sample Number');ylabel('Rect Genliği');plot(filter_y,'r');

RESULT

-----------------------------------EDIT-----------------------------------

thanks to robert bristow-johnson, i tried zero-padding to myData , end and begin.
This is i try

MATLAB CODE

% % EMG Soleus
clc;clear all;
emg_healthy=load('emg_healthy.txt');
x= 10;
fs= 4000;
my_counter = 1000;
myRawData = emg_healthy(:,2);
myRawData = myRawData(1:my_counter);
myDataRecti = abs(myRawData);
[b,a]= butter(6,x/(fs/2),'low');
for i = 1:1:(my_counter/4)
myZeroDataRecti(i) = 0;
end
for i = (my_counter/4+1):1:(my_counter+my_counter/4)
myZeroDataRecti(i) = myDataRecti(i-my_counter/4);
end
for i=(my_counter+my_counter/4+1):1:(my_counter+my_counter/4+my_counter/4)
    myZeroDataRecti(i) = 0;
end
filter_y=filtfilt(b,a,myZeroDataRecti);
firstFilter = filter(b,a,myZeroDataRecti);
flip_firstFilter = fliplr(firstFilter);
secondFilter = filter(b,a,flip_firstFilter);
flip_secondFilter = fliplr(secondFilter);
xlabel('Sample Number');ylabel('Amplitude');plot(myZeroDataRecti,'g');
hold on ;
xlabel('Sample Number');ylabel('Amplitude');plot(flip_secondFilter,'b');
xlabel('Sample Number');ylabel('Amplitude');plot(filter_y,'r');
%

RESULT

another i try , i again add zero-padding begin and end but i used filtfit method to apply without zero-padding data

MATLAB CODE

% % EMG Soleus
clc;clear all;
emg_healthy=load('emg_healthy.txt');
x= 10;
fs= 4000;
my_counter = 1000;
myRawData = emg_healthy(:,2);
myRawData = myRawData(1:my_counter);
myDataRecti = abs(myRawData);

[b,a]= butter(6,x/(fs/2),'low');
filter_y=filtfilt(b,a,myDataRecti);
for i = 1:1:(my_counter/4)
myZeroDataRecti(i) = 0;
end
for i = (my_counter/4+1):1:(my_counter+my_counter/4)
myZeroDataRecti(i) = myDataRecti(i-my_counter/4);
end
for i=(my_counter+my_counter/4+1):1:(my_counter+my_counter/4+my_counter/4)
    myZeroDataRecti(i) = 0;
end
firstFilter = filter(b,a,myZeroDataRecti);
flip_firstFilter = fliplr(firstFilter);
secondFilter = filter(b,a,flip_firstFilter);
flip_secondFilter = fliplr(secondFilter);
xlabel('Sample Number');ylabel('Amplitude');plot(flip_secondFilter,'b');
hold on ;
xlabel('Sample Number');ylabel('Amplitude');plot(myDataRecti,'g');
xlabel('Sample Number');ylabel('Amplitude');plot(myZeroDataRecti,'c');
xlabel('Sample Number');ylabel('Amplitude');plot(filter_y,'r');
%

RESULT

in my opion, there are two things we have to think about that,

1.First picture in the edit part,my signal (myCode-LEGEND) amortize itself , but matlab function (filtfilt) keep continiues , why is that ?

2.Second picture in the edit part, my signal(myCode-LEGEND) different from filtfilt function signal. And between myCode and 250Zero+originalData+250Zero i think it doesnt any phase shift. If we say that it's true what should i do ? Will i remove 250 zero from begin and end ?

As far as I can see you're doing the basic filtering right. The difference with Matlab's filtfilt function is that the latter uses certain initial conditions in the calls to filter to minimize transients.

The Matlab documentation says the following:

  filtfilt minimizes start-up and ending transients by matching initial conditions, and you can use it for both real and complex inputs

The appropriate reference is

  Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992