FFT one-sided and Parseval theorem I'm trying to get Parseval theorem working on a one-sided FFT. So far I have this code (matlab):

%% ODD
a = [1 2 3 4 5];
A = fft(a);
A1 = abs(A);
A2 = abs([A(1) 2*A(2:3)]);

Ea=sum(a.^2)
EA1=sum(A1.^2)/5
EA2=sum(A2.^2)/5

%% EVEN
b  = [1 2 3 4 5 6];
B  = fft(b);
B1 = abs(B);
B2 = abs([B(1) 2*B(2:3) B(4)]);

Eb =sum(b.^2)
EB1=sum(B1.^2)/6
EB2=sum(B2.^2)/6

And the output is:

Ea =  55
EA1 =  55
EA2 =  65
Eb =  91
EB1 =  91
EB2 =  107

The scaling for the one-sided is to preserve amplitude information and not for calculating Parseval theorem?

Parseval's theorem says that the the following relationship holds
$$
\sum_{n=1}^{N} a[n]\,a^*[n] = \frac{1}{N}\sum_{k=1}^{N} A[k]\,A^*[k]
$$
where $A[k]$ is the discrete Fourier transform of $a[n]$, both assumed to be of length $N$ (no padding). This arises from the fact that the signal energy calculated from the time domain and frequency domain must be equal. See this answer for a little more detail about how this formula comes about.
If you want to calculate a one-sided spectrum, you need to throw away the negative frequencies (which are redundant when the signal is real). However, because half of the energy was contained in these entries, you need to multiply the energy in the remaining bins by 2 - except for the DC and Nyquist. I see that you have recognised correctly that the Nyquist is only present when $N$ is even, and that the first entry is the DC component regardless of $N$ - good. However, you have multipied the amplitudes by 2 (straight after the FFT), not the energy - the negative frequency bins contain half of the energy.
Instead, once you have calculated the discrete Fourier transform using the FFT algorithm (your variables $A$ and $B$), I would first then obtain the energy spectral density of the signal $a[n]$ using
$$
\textrm{ESD}_a[k] = |A[k]|^2 = A[k]\,A^*[k]
$$
and then throw away the negative frequencies from this vector instead. This way you throw away half of the energy correctly, and you can confidently multiply the energy in the remaining bins by $2$.
The following code illustrates this:
%% ODD
a = [1 2 3 4 5];
A = fft(a);
ESD_a = A.*conj(A);
ESD_a_onesided = [ESD_a(1) 2*ESD_a(2:3)];

E_a_timedomain = sum(a.^2)
E_a_twosided = sum(ESD_a)/5
E_a_onesided = sum(ESD_a_onesided)/5

%% EVEN
b = [1 2 3 4 5 6];
B = fft(b);
ESD_b = B.*conj(B);
ESD_b_onesided = [ESD_b(1) 2*ESD_b(2:3) ESD_b(4)];

E_b_timedomain = sum(b.^2)
E_b_twosided = sum(ESD_b)/6
E_b_onesided = sum(ESD_b_onesided)/6

The result is then correctly
E_a_timedomain = 55
E_a_twosided = 55
E_a_onesided = 55

E_b_timedomain = 91
E_b_twosided = 91
E_b_onesided = 91

EDIT----------------------------
Actually, the energy values of $55$ and $91$ obtained are only correct if we are assuming the sampling period of the signal acquisition was $T_s=1$.
The signal energy of a continuous signal $a(t)$ is defined as
$$
E_s = \int_{-\infty}^{+\infty}|a(t)|^2\;dt
$$
and the energy of a sampled version of it is then
$$
E_s = \sum_{n=1}^{N}|a[n]|^2\cdot T_s
$$
and you can see that we need to account for the signal sampling period to get the right energy.
The following verifies Parseval's theorem if we for example had collected the signal at some other sampling period $T_s\neq 1$
Ts = 0.05; % the sampling period of acquisition
Fs = 1/Ts; % the sampling frequency of the acquisition

%% ODD
a = [1 2 3 4 5];
N = 5

A = fft(a)*Ts;
ESD_a = A.*conj(A);
ESD_a_onesided = [ESD_a(1) 2*ESD_a(2:3)];

E_a_timedomain = sum(a.*conj(a))*Ts
E_a_twosided = sum(ESD_a)*Fs/N
E_a_onesided = sum(ESD_a_onesided)*Fs/N

%% EVEN
b = [1 2 3 4 5 6];
N = 6

B = fft(b)*Ts;
ESD_b = B.*conj(B);
ESD_b_onesided = [ESD_b(1) 2*ESD_b(2:3) ESD_b(4)];

E_b_timedomain = sum(b.^2)*Ts
E_b_twosided = sum(ESD_b)*Fs/N
E_b_onesided = sum(ESD_b_onesided)*Fs/N

with the output
E_a_timedomain = 2.75 [signal^2 sec]
E_a_twosided   = 2.75 [signal^2 sec]
E_a_onesided   = 2.75 [signal^2 sec]

E_b_timedomain = 4.55 [signal^2 sec]
E_b_twosided   = 4.55 [signal^2 sec]
E_b_onesided   = 4.55 [signal^2 sec]


changing filtering "speed" (cutoff frequency) of an IIR filter knowing only its coefficients i got an IIR filter bu I have only the coefficients.
now i'd like to be able to change the filter characteristics of the low pass filter (the "cutoff" frequency) but all i got are a and b coefficients.
i'd like  to be able to set a multiplier that will "scale" the filter by that amount. if i set the multiplier to 0.5 i'd like to filter the same but two times lower or two times slower (50 Hz cutoff instead of 100Hz)
is it possible ? Or is it a lost cause ?
and with Matlab ?
thanks so much
Jeff

Your probably your best option here is "frequency warping". You can calculate the poles and zeros from the coefficients, warp the poles and zeros and then recalculate the coefficients again.
Warping is a procedure the applies a conformal mapping to the poles/zeros. Specifically a conformal mapping that maps the unit circle onto itself will maintain the overall shape of the filter but either bunch it up or stretch in the low frequencies (and vice versa in the high frequencies). A good choice is a first order allpass filter, that indeed maps the unit circle onto itself.
Loosely speaking you replace all delays in your filter with a first order all pass filter and recalculate the poles and zeros. That takes a bit of math to work out the details, but the actual code would be quite efficient.


What factors go into choosing DCT type 1 over type 2 or type 3 etc? I'm a newcomer to the DCT, DFT etc. I understand that they're intrinsically related, I understand that the DCT tends to "compress" the information in a signal towards the first few bins, leaving the rest relatively close to zero. I'm trying to learn more.

DCT has different "types". Based on reading wikipedia, the types mainly differ in adding 0.5 to either the input index, the output index, or both, when computing the matrix coefficients. Type 1 does neither, type 4 does both, and 2/3 do one or the other.

It's not clear why types 2, 3, or 4 would ever valuable. What benefit does shifting the input or output signal by half give? What are some examples of applications that are well-suited for DCT type 2,3 or 4, but are not well suited for DCT type 1? If I have an array of floats that I want to compute a DCT on, what factors should go into an educated decision on which type to use?

As you described correctly, DCT is closely related to DFT. The two main motivations behind DCT are it's being a real transform and its compaction of signal energy. The first makes it computationally simpler, while the second makes it a prime tool for transform based lossy data compression applications. As you understood correctly.

Then what remains is where the types come from and why is it the DCT-II as the most common choice for data compression? The answer is as follows:

Given a signal $x[n]$ of length $N$, The different DCTs are defined by forming even-symmetric extensions of shifted copies $x[n]$ in a number of possible ways, each different combination yielding a different type of DCT. Also the odd-symmetric extensions are called as the DST (Discrete Sine Transform).

The DCT-II , by virtue of its type of symmetry, has the most energy compaction among other DCT types and therefore is naturally the choice.

The DCT-II signal formation (based on A.Opp DTSP-2e) is as follows:

Consider a signal $x[n]$ of length N, whose DCT-II we want to compute. First extend $x[n]$ by padding $N$ zeros and create the extended signal $x_e[n]$

Then the DCT Type-II signal $\tilde{x}_2[n]$ is a periodic signal with period $2N$ and related to $x_e[n]$ as:
$$\tilde{x}_2[n] = x_e[((n))_{2N}] + x_e[((-n-1))_{2N}]$$
whose base period $x_2[n]$ extends from $n=0$ to $n=2N-1$. And where the $((\cdot))_{N}$ operation here is the modulus operator, providing the periodicity of $\tilde{x}_2[n]$.

Eventually, it is shown that the $N$-Point  unnormalised  DCT-II, $X_{c2}[k]$, of the signal $x[n]$ is related to the $2N$-Point DFT, $X_2[k]$, of the signal $x_2[n]$ as:
$$X_{c2}[k] = e^{-j\pi k/{2N}} X_2[k]$$
for $k=0$ to $k=N-1$

The following Matlab script creates the Type-II DCT signal $x_2[n]$ from $x[n]$ and finds DCT of x[n] from DFT of x2[n] for comparison

N=8;                            % Period for x[n]
n = [0:2*N-1];                  % A set of sequence indices of at least length 2N
x = [1:N];                      % A ramping test sequence for x[n] of period N
xe = [x zeros(1,N)];            % Extended sequence from x[n], by padding N zeros
x2 = xe(1 + mod(n,2*N)) + xe(1 + mod(-n-1,2*N));    % DCT-II type periodic signal
DFTx2 = fft(x2,2*N);            % 2N-Point DFT of x2[n]
Dft2Dct = real( exp(-j*pi*[0:N-1]./(2*N)).*DFTx2(1:N) );   % Create N-DCT of x[n] from 2N-DFT of x2[n]
DCTx = sqrt(2*N).*dct(x,N);     % N-POint DCT of x[n], for test purposes
DCTx(1)= sqrt(2)*DCTx(1);       % UNnormalise the Matlab's default normalised DCT-II


Image processing - Is mean filter shift-invariant? Update:
So, according to this article.
For a 2D mean filter,
$$
H[f(x,y)]=\frac{1}{NM}\sum_{k=0}^{N-1}\sum_{p=0}^{M-1}f(x-k,y-p)=g(x,y)\tag{1}
$$
$$
g(x+x_0,y+y_0)=\frac{1}{NM}\sum_{k=0}^{N-1}\sum_{p=0}^{M-1}f(x+x_0-k,y+y_0-p)\tag{2}
$$
$$
H[f(x+x_0,y+y_0)]=\frac{1}{NM}\sum_{k=0}^{N-1}\sum_{p=0}^{M-1}f(x+x_0-k,y+y_0-p)\tag{3}
$$
because (2) equals to (3), a 2D mean filter is shift-invariant.
Is the proof right?
BTW, if f(x,y) is a 2D image, such as

1 2 1
2 3 2
1 2 1

and g(x,y) is (using zero padding)

 8/9 11/9  8/9
11/9 15/9 11/9
 8/9 11/9  8/9

What does these
$$
H[f(x+x_0,y+y_0)]\tag{4}
$$
and
$$
g(x+x_0,y+y_0)\tag{5}
$$
mean?
(I want to prove it by using a image.)
Origin:
I read this article
and this Wikipedia page.
However, I still don't know what it means for (or how to apply to) a mean filter.
Can anyone give me an image example?
(for example, https://dsp.stackexchange.com/a/14435 is easier to understand)

The 1D mean filter defined by the following I/O relation is time/shift invariant;

$$y[n] = T\{x[n]\} = \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d_1}^{n+d_2}{x[k]}$$

For $d_1$ and $d_2$ being positive integers, the filter computes the mean (or the average) of its input $x[n]$ in the range of its summation index. To show that it is shift invariant we shall simply show the truth of the following statement:  $y[n-d] = T\{x[n-d]\}$ for and $d$ and any input $x[n]$

$$ y[n-d] = \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d-d_1}^{n-d+d_2}{x[k]} ~~~\overset{?}{=}~~~ T\{x[n-d]\} = \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d_1}^{n+d_2}{x[k-d]}$$

Which is easily shown to be true by making the simple change of variable in the second sum as $k-d=k'$ and then optionally replacing $k'$ with $k$ again:
\begin{align}
T\{x[n-d]\} &= \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d_1}^{n+d_2}{x[k-d]} ~~~~ &\scriptstyle{\text{replace k-d with k'}}\\
T\{x[n-d]\} &= \frac{1}{d_1 + d_2 + 1} \sum_{k'=n-d-d_1}^{n-d+d_2}{x[k']} &\scriptstyle{\text{replace k' with k again}}\\
T\{x[n-d]\} &= \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d-d_1}^{n-d+d_2}{x[k]} &\scriptstyle{\text{recognise that this is =  y[n-d]}}\\
T\{x[n-d]\} &= y[n-d] &\scriptstyle{\text{Hence the given System is LTI}}\\
\end{align}

This is extended to 2D in a straight-forward manner


What is the largest "safe" order for the digital Butterworth filter of a given signal? Say I want to apply a zero-phase Butterworth bandpass filter, implemented in MATLAB using butter and filtfilt. So, for an order-$4n$ filter with passband $[f_{\mathrm{min}\,},f_{\mathrm{max}}]$, I believe the transfer function $H$ is given by
$$ H(2\pi if) \ = \ \frac{1}{1+\left( \frac{f}{f_\mathrm{max}-f_\mathrm{min}} + \frac{f^{-1}}{f_\mathrm{max}^{-1}-f_\mathrm{min}^{-1}} \right)^{\!2n}}\,. $$

    It seems that for any given digital signal $x(t)$ (other than the constant zero signal), there is a critical number $N_x$ such that:

    for any $n < N_x$ the results are fairly sensible;
    but for $n \geq N_x$, the results change considerably and are unreliable.

    (1) What is the cause of this behaviour? Is it due to numerical imprecisions that grow large when the stability of the filter is poor?

    (2) Is there a straightforward way to determine what this critical number $N_x$ is?

Let me illustrate the phenomenon with a couple of examples:

EXAMPLE 1

Consider the signal $x(t)=\sin(2\pi t)$, $t \in \{0,0.01,0.02\ldots,99.99,100\}$. (This signal has 10001 points altogether.)

Suppose I apply the Butterworth filter with passband $[f_{\mathrm{min}\,},f_{\mathrm{max}}]=[0.5,3]$. Of course, I expect my filtered signal $x_{\mathrm{filt}}(t)$ to be basically the same as $x(t)$, except:

there may be edge effects near the start and/or end;
for $n=1$ or $2$, the amplitude is likely to be slightly reduced, as the frequency of the sinusoid $x(t)$ is not very close to the centre frequency $\sqrt{1.5}$ of the filter.

For each $n \leq 5$, the results are as expected. E.g., for $n=5$ (corresponding to filter order 20), the result is as follows:

However, when I now take $n=6$ (filter order 24), I suddenly get a drastically different result:

A wavelet transform of the original and filtered signal for $n=6$ looks as follows:

For even larger $n$, the results simply blow up. For $n=7$, the filtered signal is equal to NaN for about the first two thirds of the duration, and then starts taking values of the order of $\,-10^{305}$, going up to a maximum value of about $4 \times 10^{231}$.

So in this example, $N_x=6$.

EXAMPLE 2

This is exactly the same as Example 1, except with the sampling frequency considerably reduced.

Consider the signal $x(t)=\sin(2\pi t)$, $t \in \{0,0.08,0.16\ldots,99.92,100\}$. (This signal has 1251 points altogether.)

Once again, suppose I apply the Butterworth filter with passband $[f_{\mathrm{min}\,},f_{\mathrm{max}}]=[0.5,3]$. For $n \leq 16$, I get sensible results. Here is the filtered signal with $n=16$:

I think the "early warning sign" here is that the slight inaccuracy of the filtered signal persists throughout the whole duration, not just at the start and end; for $n \leq 15$, the visible inaccuracy is only at the edges.

For $n=17$, I get quite a bad result:

In this case, Fourier and wavelet transforms are as follows:

Once again, for $n=18$, the results blow up:

(Please note that the values of the filtered signal are extremely large throughout, not just for $t<35$. At $t=70$, the value is about $-6.5 \times 10^{11}$.)

So in this example, $N_x=17$.

REMARK: In many (probably most) other examples, there is no intermediary $n$-value between the $n$-values giving consistent results and the $n$-values for which the filtered signal blows up to extreme values. In this regard (but only this regard), the above two examples seem somewhat exceptional.

One cause is that higher order Butterworth filters have poles closer to the unit circle.  This nearby infinite gain point increases the likelihood of numerical instabilities.  (e.g. rounding/arithmetic/quantization noise may move a pole to the “wrong” side of the unit circle.)
Where the numerical noise will blow up depends on your executable code’s exact arithmetic sequence, and the quantization and number format used.  Often, dividing up a high-order filter into a sequence of biquad stages produces a more stable (better conditioned) sequence of arithmetic operations.


Perfect reconstruction in two-channel filter bank [from Discrete-time Signal Processing by Oppenheim and Schafer, 3rd ed.]
I was not able to prove this:

If the analysis and synthesis filters are ideal so that they exactly split the band 0 ≤ |ω| ≤ π
into two equal segments without overlapping, then it is straightforward to verify that
Y(ejω) = X(ejω); i.e., the synthesis filter bank reconstructs the input signal exactly.

$G_0(e^{j\omega})$ and $H_0(e^{j\omega})$ are ideal low pass filters with passband $[0,\pi/2]$ and stopband $[\pi/2,\pi]$. $G_1(e^{j\omega})$ and $H_1(e^{j\omega})$ are ideal high pass filters that are complementary to the lowpass filters, i.e., their passbands coincide with the lowpass filters' stopbands and vice versa.
Note that shifting the frequency axis by $\pi$ transforms a lowpass filter to a highpass filter and vice versa. I.e., $H_0(e^{j(\omega-\pi)})$ is a highpass filter, and $H_1(e^{j(\omega-\pi)})$ is a lowpass filter. Consequently, $G_0(e^{j\omega})H_0(e^{j(\omega-\pi)})=0$ and $G_1(e^{j\omega})H_1(e^{j(\omega-\pi)})=0$. So the second term $(111b)$ corresponding to the aliased input spectrum vanishes. We're left with the first term, which is constant because we add an ideal lowpass spectrum and its complementary highpass spectrum.


Latency handling at filter output In MATLAB, latency at the start of the output of a filter can be handled by discarding those bits and starting the signal where the latency ends. but how can the garbage value at the end be handled?

You may use the internal state variable zf that is available as an output from the filter function

[y,zf]=filter([ 1 2 3],1,[1 0 0 0 1 ])

it contains the final condition of the filter, in the example above (with b=[1 2 3] and a=1) you will get

y =

     1     2     3     0     1

zf =

     2
     3

Concatenating zf to the end of y will give you the complete output signal, it is the same as adding zeros to the input signal at the end.

 y=filter([ 1 2 3],1,[1 0 0 0 1 0 0])

 y =
     1     2     3     0     1     2     3


non-uniform sampling on the unit circle - applications If there is an N-point sequence, it has a Fourier Transform.  DFT of this sequence is the sampling of the Fourier Transform at N equally spaced points on the unit circle.  What happens if we take some M non-uniformly spaced samples on the unit circle with M < N, M = N and M > N.  Any studies done on this.  If yes, what are the applications?

Uniform sampling is the definite standard as it simplifies not only mathematical derivations but also practical applications. On the other hand, there are some applications (such as array of sensors located in a nonuniform grid) in which taking uniformly distributed samples from a signal of interest is impossible. In such cases nonuniform sampling in either time or space is performed and the usual approach is to convert those nonuniform samples to uniform ones. Then follows the standard operations on those precomputed uniform samples.


How to perform median filtering in multi-label prediction I'm trying to perform automatic pitch tracking for piano. It's easy to detect the note when there's only a single note pressed. There are 88 keys on the piano and I set the result to one of them. There may be misclassifications due to octave shift. For example, instead of 40, I find 52. They are both C, but 52 is one octave up. I can remove most of these octave shifts by performing median filtering on the results.

But it is confusing when I move to polyphonic pitch tracking. When there are more than a single note pressed together, there can be varying number of pitches detected at time t.

For example, let labels(f,t) be 1 if note f is detected at time t and let the binary matrix labels it be as follows

labels = 0 0 0 0 1 ...
         1 1 0 1 0 ...
         0 1 0 0 1 ...
         .............
         0 0 1 0 1 ...
         .............
         0 0 0 0 1 ...

Let the detected notes be as follows correspondingly:

labelsAsNotes = 2; 2 3; 14; 2; 1 3 14 88; ...

How can I do median filtering in this case? Do you have any other recommendations to remove octave shifts?

I don't think median filtering is a good approach for this problem. It is definitely used as a smoothing mechanism in single pitch tracking, but I don't think that's necessarily the best way to go. For example, when tracking single pitches, a dynamic programming algorithm can be really powerful. At each frame of audio you have pitch candidates with some strength to them. Dynamic programming lets you find an optimal path between those candidates such that you penalize big jumps in pitch (like a whole octave) while promoting pitch candidates with high weights.

As far as I'm aware, though, this doesn't scale well to multiple pitches. I have seen authors use this algorithm to track just two pitches in the case of two concurrent speakers.

This is a really hard problem to solve in general, and I don't think just filtering that binary matrix will be sufficient -- I think you want some sort of weights associated with each pitch candidate, and to make decisions with those weights. You can also take time into account -- if you detect two notes that are an octave (or other integer multiple) apart, and follow the same relative pitch trajectory and onset at exactly the same time, they are probably the same note. You could also try using the weights to inhibit octave errors -- for example, if you detect a note, but there's significant energy at half its frequency, you might be looking at an octave error.


Is Linear Prediction ever exact? I'm just learning about linear prediction and was wondering whether or not there are any signals that would produce no error if run through a linear prediction algorithm. I was thinking (maybe naively) that if we have a linear signal (one that changes linearly with time) that we should be able to linearly predict this with no error. Is my intuition on this correct? Are there any other signals that are like this?

First of all, thanks so much for the replies! Secondly, I'm glad my intuition was correct, but I was having doubts in the first place because I couldn't exactly predict future samples of my signal no matter what my input signal was. For example, here's an image of me trying to get the linear prediction coefficients using 'lpc' function in Matlab.

I've tried this for several of the types of signals suggested above, and didn't get the correct coefficients for any of them. Why does this happen? I doubt the 'lpc' function is incorrect so it just must be doing something different than what I thought it was doing.

Theoretically, yes. So you can for instance predict values taken by polynomials, with the following initial patterns:

degree $0$: $\hat{S}[n] = {S}[n-1]$;
degree $1$: $\hat{S}[n] = 2{S}[n-1]-{S}[n-2]$;
degree $2$: $\hat{S}[n] = 3{S}[n-1]-3{S}[n-2]+{S}[n-3]$;
degree $3$: $\hat{S}[n] = 4{S}[n-1]-6{S}[n-2]+4{S}[n-3]-{S}[n-4]$;

and so on. You can recognize binomial coefficients.

There are similar  linear relations for other signals when sampling conditions are met, and for different cost functions. I am happy to dig into old stuff here. In a 1927, G.  Udny Yule (the same as in autoregressive modeling known under Yule-Walker), On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers,  gave the following trigonometric recurrence relation  (using the same notations), for:

$$\begin{align}
u_0 & =& A \sin\left({2\pi\frac{t+\tau}{T}}\right)\\
u_1 & =& A \sin\left({2\pi\frac{t+\tau+h}{T}}\right)\\
u_2 & =& A \sin\left({2\pi\frac{t+\tau+2h}{T}}\right)\end{align}
 $$
then
$$ u_2 = (2-\mu)u_1-u_0$$
with $\mu = 4 \sin^2 \pi\frac{h}{T}$. One can notice that, when $\mu =0$, one recognizes the linear extrapolation: $ u_2 = 2u_1-u_0$.

I would suggest:

The Theory of Linear Prediction by P. P. Vaidyanathan, or
chapter 13.6 Linear Prediction and Linear Predictive Coding  in Numerical Recipes.

In practice, the exactness is limited by the data precision (in time and amplitude), for instance by finite-word quantization (rounding errors, floating point operations) and jitter, in time.


Is information lost due to Hanning window perturbations? I apologize for the obscure phrasing of the title of the post...My question is quite simple...

Let us imagine an audio signal (any song will do). We want to apply the Short Time Fourier Transform algorithm, so we break our signal into frames and view each frame through a window. For the sake of the example let us assume the window is a Hanning window.

Many algorithms have a parameter to specify the sample at which to center the Hanning curve in a given frame. Now, let us assume that the dominant signal in that frame happens to fall near one of the the tails of the Hanning window for that frame. The window will then squash the importance of the sample and boost some samples which were not that prominent in the original signal. In essense, we have lost some possibly valuable information.

Would I then be correct in assuming that not only is the number of frames of importance in FT, but also the samples at which the window functions are centered?

Thanks in advance.

The STFT itself is invertible, which means by means of the STFT no information about the original signal is lost.

For common implementations of the discrete STFT, we do not calculate the FFT of each possible shift of the window (i.e. every sample), but instead we have a bigger step size (of a few samples in time domain). This sampled STFT (because we evaluate the STFT(t,f) only at t=k*T, where T is the spacing between the windows) is actually a discrete Gabor Transform (of a 1D signal, most references online refer to Gabor Transform of images). In superficial terms, the Gabor transform is invertible as long as the sampling of the STFT is not too coarse. A necessary condition for invertablity (in continuous time) is that $TF\leq 1$, where $T$ is the sampling distance in time and $F$ is the sampling distance in frequency (since in continuous-time, we have infinite frequency resolution).  The invertability becomes the easier, the smaller $TF$, i.e. the finer the sampling of the STFT is done.

So, to answer your question: Despite sounding intuitive that some portion of the signal is lost, if the window is positioned at "unlucky" positions, the STFT does not lose information about the original signal, as long as it is sampled fine enough. So, the signal itself does not determine, if information is lost, but the sampling of the STFT does.


Bessel filter second-order sections Q and Fc multiplier derivation A practical guide for digital IIR audio filters has cookbook-style values for creating higher-order Bessel filters out of biquads, but the values listed aren't very precise:

You should multiply the Fc for each stage by the following coefficients:
1: 1.00
2: 1.27
3: 1.32 1.45
4: 1.60 1.43
5: 1.50 1.76 1.56
6: 1.90 1.69 1.60

The corresponding Q values are:
1: ---
2: 0.58
3: --- 0.69
4: 0.81 0.52
5: ---- 0.92 0.56
6: 1.02 0.61 0.51

What are the expressions to calculate these values exactly?

For the Butterworth table, the values are given by Q = 1/(2 * sin((pi / N) * (n + 1 / 2))), for instance.  So where the table says 4: 0.54 1.31, the equation gives 1.3065629648763766 and 0.54119610014619701.  But I can't find the expression for Bessel filters, except that 2: 0.58 is $1\over\sqrt3$.

Prototype-generating functions just have long tables of numbers, so maybe these are not easy to calculate?  (Bessel functions?)  If so, a listing of these Q and Fc values for the first several orders would be good.

(And yes, I know the bilinear transformed Bessel is no longer linear-phase.  :/  I just want to know where these numbers are from.)

There is no simple closed-form equation to calculate Bessel filter bi-quad coefficients. The poles of the filter come from a Bessel polynomial. Higher order Bessel polynomials are determined using a recursion relationship. You need computer algorithms to calculate the poles and resolve the bi-quad coefficients.

Two books I can recommend: “Passive and Active Filters” by Wai-Kai Chen and “Analog Filter Design” by M.E. Van Valkenburg.


Question about aliasing As far as I understand, you can have two different continuous-time signals with the same discrete-time frequency spectrum after they are sampled and it may be possible these shifted replicas in the discrete frequency domain do not overlap with each other. In textbooks aliasing is associated with overlapping of shifted replicas. is the case I described a case of aliasing?

What you describe is indeed aliasing. If there are frequency components above half the sampling frequency then they will be folded back to the interval $[0,f_s/2]$. This folding back is called aliasing, no matter whether these aliased components overlap with other parts of the signal spectrum or not.
Note that you can also use aliasing to your advantage if you're careful that no overlap occurs. An example is bandpass sampling, where you deliberately choose a sampling frequency that is much smaller than twice the highest signal frequency.


Don't understand the relation about the wireless communication flow chart and TDMA illustration diagram The first picture is the wireless communication flow chart,and the second is the TDMA illustration diagram,Which part of first picture is TDMA ? or if first picture is the standard process for TDMA,then how to explain it by using the TDMA illustration diagram? I mean, TDMA illustration diagram just show me  that divide whole time into several part ,and it seems that there is no relation about the block in the first picture like source coding,modulation ,etc.

TDMA wraps around the outside of your block diagram of the communication link.  There are $M$ users who are buffering the data that they wish to transmit. The TDMA switch at the input selects each user in turn and connects its buffer to the block marked Source Coding for a fixed amount uf time. During this time, the contents of the user's buffer are source coded, channel coded, transmitted etc.  When the user's time is up, the TDMA switch connects the next user's buffer to the Source Coding block, etc.  At the receiver end, there is a synchronous switch that routes the Source Decoder's output to the appropriate user's buffer which the user empties at a lower rate (if it so chooses so that it gets its data at a constant rate rather than in bursts).


Sampling Frequency Math The Nyquist formula for digitally sampling and recovering a continuous waveform is well known. Are there any mathematical formulas that show that it is computationally infeasible (or impossible) to reconstruct a waveform if you go below a certain sampling frequency?

Sampling a continuous signal $x(t)$ can be modeled by multiplication with a delta pulse train:

$$x_s(t)=x(t)\cdot\sum_{n=-\infty}^{\infty}\delta(t-nT)=
\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT)$$

where $T$ is the sampling interval: $T=1/f_s$ ($f_s$ is the sampling frequency). Multiplication in the time-domain corresponds to convolution in the frequency domain. So the spectrum of the sampled signal can be computed by convolving the transforms of $x(t)$ and of the pulse train:

$$x(t)\Longleftrightarrow X(f)\\
\sum_{n=-\infty}^{\infty}\delta(t-nT)\Longleftrightarrow
\frac{1}{T}\sum_{k=-\infty}^{\infty}\delta(f-kf_s)
$$

Convolution of these two Fourier transforms gives

$$x_s(t)\Longleftrightarrow X_s(f)=
\frac{1}{T}\sum_{k=-\infty}^{\infty}X(f-kf_s)\tag{1}$$

Equation (1) shows that the spectrum of a sampled signal is composed of shifted versions of the original spectrum. The shift corresponds to multiples of the sampling frequency $f_s$. Consequently, if the spectrum $X(f)$ is exactly zero for frequencies greater than $f_s/2$, then the shifted versions do not overlap, and the original spectrum can be reconstructed by lowpass filtering the spectrum of the sampled signal. If $X(f)$ is not band-limited to $f_s/2$, i.e. it has frequency components greater than $f_s/2$, the shifted spectra in Equation (1) overlap, and the original signal can in general not be reconstructed from the sampled signal. This is what the sampling theorem states.

EDIT: For non-periodic signals - i.e. if $X(f)$ does not contain Dirac impulses - it is sufficient for reconstruction that $X(f)$ does not contain frequency components above $f_s/2$. In this case it is (theoretically) possible to reconstruct $x(t)$ from its samples $x(nT)$ by applying an ideal lowpass filter with cut-off frequency $f_s/2$. The impulse response of such an ideal lowpass filter is given by

$$h_{LP}(t)=\frac{\sin (\pi t/T)}{\pi t/T}$$

Filtering $x_s(t)$ with this impulse response yields the original signal

$$x(t)=\sum_{n=-\infty}^{\infty}x(nT)\frac{\sin [\pi (t-nT)/T]}{\pi (t-nT)/T}\tag{2}$$

Equation (2) is usually referred to as the Whittaker-Shannon interpolation formula.

Normally one does not need to consider what is actually happening at the frequency $f_s/2$. The only case where this is necessary is the presence of a pure sinusoid with frequency $f_s/2$. From Equation (1) we can see that in this case the shifted spectra (for $k=\pm 1$) will also produce spectral lines at $\pm f_s/2$, which in general will make impossible the reconstruction of the sinusoid at $f_s/2$. Apart from this somewhat artificial case, it suffices to require that the signal does not contain any frequency components greater than $f_s/2$.

The above considerations are of course of a purely theoretical nature, because there exist neither ideally band-limited signals nor ideal lowpass filters. Consequently, in practice we will choose the sampling frequency such that it is well above twice the maximum signal frequency of interest. The anti-aliasing filter must pass frequencies up to this maximum frequency, then it has a finite transition bandwidth, and then we can choose a frequency in the stopband of the filter where the stopband attenuation guarantees that aliasing errors are sufficiently small for the given application. This frequency can finally be considered the "bandwidth" $B$, and the sampling frequency can now be chosen as $f_s\ge 2B$.


If a filter has a non-linear phase lag will the output signal be distorted? EDIT: From here (https://www.mathworks.com/help/signal/ref/butter.html)
it appears that an Butterworth filter has a non-linear phase lag at the normalized frequency of 0.5 onwards.
Does that mean that if the signal has multiple frequencies beyond 0.5, its output after passing through this filter will be distorted by these frequencies?

NOTE: To help clarify things, this question had been modified from the original one below.

A butterworth low pass filter would have a large phase lag for signals passing through at higher frequencies.
[EDIT:
The assumption is, for some filters to work, you'd have to do some kind of transform on the input (into the s domain), multiply that by the filter's s-domain transfer function (such as ones from here ece.uic.edu/~jmorisak/blpf.html), then to see the time-domain filtered output, the resulting output has to undergo an inverse transform.]
I envision that this means that the output signal coming out of the filter, when converted back into time domain, would come out distorted due to this lag.
[EDIT: The assumption is, phase lag is the same as a time delay between the input and the output, and if the input has more than one frequency component, the different components would be delayed differently, and end up being combined with another signal output at a later point in time.]
Would there be similar distortions (i.e.: as if there were a phase lag) for an IIR time-domain digital filter?

I think your question comes from several misunderstandings. The fact that the phase lag of a system becomes more negative for large frequencies does not mean that there's more distortion of larger frequencies. Neither does it mean that high frequency signal components experience more delay when passing through the filter.

Imagine an ideal system that simply delays the input signal while keeping its amplitude unchanged. So for an input signal $x(t)=\sin(\omega t)$ the output signal is given by

$$y(t)=\sin(\omega(t-\tau))$$

where $\tau$ is the delay. The output signal $y(t)$ can be written as

$$y(t)=\sin(\omega t-\omega\tau)=\sin(\omega t+\phi(\omega))$$

where $\phi(\omega)=-\omega\tau$ is the phase response of the system. Note that the phase becomes more negative for larger frequencies, and this is exactly what is needed for a system that does not distort the input signal.

What causes distortion are deviations of the system's magnitude response from the value $1$, and deviations of the system's phase response from the ideal linear phase shown above.

And, as mentioned in the comments, infinite impulse response (IIR) filters like a Butterworth filter are almost always implemented in the time domain. You implement a given difference equation, which is equivalent to implementing a system with a given impulse response and a given frequency response (the latter two are Fourier transformation pairs), and the phase response is part of the system's frequency response.


Sinusoidal Wave Terminology Suppose I have a sinusoidal signal of the form $x(t) = a\sin\left(2\pi b(t-c)\right) + d$.  I understand $a$ is the amplitude, $b$ is the wave number, and $c$ is the phase shift.  I know that the term $d$ shifts the entire signal vertically, but does it have a specific name in signal processing?

It is usually called DC offset or DC component. Since you use $t$ as a variable, which suggests "time", the quantity $b$ would normally be called frequency. Also note that $c$ is not the phase but a time offset. The phase $\phi$ is usually defined as in

$$x(t)=A\sin(2\pi ft+\phi)\tag{1}$$

Comparing (1) to your formula gives $\phi=-2\pi bc$.


Problem converting from Eb/N0 to SNR I'm doubtful about which of these is correct:

$\frac{E_b}{N_0}=\frac{SNR}{k}\frac{T_{symbol}}{2T_{sampling}}=\frac{SNR}{k}\frac{L}{2}$
$\frac{E_b}{N_0}=\frac{SNR}{k}\frac{T_{symbol}}{T_{sampling}}=\frac{SNR}{k}L$

$L$ is the interpolation factor and $k=log_2(M)$ (number of bits per symbol).

There's just a difference of a $1/2$ factor. I know where this relation comes from, and that $1/2$ depends on noise power being $\sigma_N^2=N_0/2$ (baseband) or $\sigma_N^2=N_0$ (passband). Since I want to generate $\sigma_N^2$ given a certain $E_b/N_o$, I'm not sure which one to pick.

If I have $y(t)=x(t)+n(t)$, with $n(t)$ being AWGN and $x(t)=s(t)\cos(\omega_c t)$, so $x(t)$ is passband and $s(t)$ is baseband, I'm not sure if I should pick the second one because I'm adding it to $x(t)$ or the first one.

I'm not comfortable with either of your two equations, in part because you don't define any of the three variables involved, but also because I think there's an easier way to understand $E_b/N_0$ and simulate it.

I think the most important thing to know about $E_b/N_0$ is that it is measured at the matched filter's output. The first consequence of this is that the measurement is done at baseband.

Let's assume that your signal $s(t)$ is linearly modulated, so $$s(t) = \sum_k a_k p(t-kT_p),$$ where $R_p = 1/T_p$ is the baud rate and the information is transmitted by the symbols $a_k$ which are elements of an $M$-PAM constellation with average energy $E_p$. Furthermore, the Nyquist pulse $p(t)$ has unit energy.

Let's assume the system has zero noise. In this case, the output of the matched filter, sampled at rate $R_p$, is equal to $a_k$ and its average energy is $E_p$. Then, the average energy per bit at the matched filter's output is $$E_b = \frac{E_p}{\log_2(M)}.$$

Now let us consider the noise. Assume that the input to the matched filter is $n(t)$, which has power spectral density $N_0/2$. Then the filter's output sampled at rate $R_p$ is a Gaussian random variable with variance $\sigma_n^2 = N_0/2$.

The signal-to-noise ratio is defined as the ratio of the average energy per bit to the average energy of a noise sample: $$\text{SNR} = \frac{E_b}{\sigma_n^2} = \frac{2E_b}{N_0}.$$ You often see the SNR defined instead as $E_b/N_0$; the reason is that, in a quadrature system, the total noise energy in the receiver is $2\sigma_n^2 = N_0$. As long as the definition is made clear, one could use either one.

Let's see how to simulate this. Say you want to simulate 4-PAM at $E_b/N_0 = 10$. First, generate 4-PAM symbols with $E_b=1$, or $E_p = 2$. Then, generate noise with variance $\sigma^2 = N_0/2 = E_b/(2 \times 10) = 1/20$. Assuming Matlab:

% generate 100 random a_k
a = 2*(randi([0,3],1,100)-1.5);
% scale them so that Ep=2, Eb=1
a = a*sqrt(2/5);
% Noise has N0 = 1/10, so N0/2 = sigma2 = 1/20
n = sqrt(1/20)*randn(1,100);
% matched filter's output:
y = a + n;

Now y has $E_b/N_0 = 10$.

Note that extending this to quadrature modulation is trivial: just think of the quadrature system as two parallel PAM systems, each with noise variance $N_0/2$.

I realize that I have only indirectly answered your question, but I hope this is still useful.


How to change x-axis from samples to time? I need to change the x-axis on this graph from 'samples' to time:

I created the graph using this process:

y, _ = librosa.load('sound/data/kea-song.mp3', 48000)
y /= y.max()

# compute the rmse
e = librosa.feature.rmse(y=y)[0]
e -= e.min()#what does this really do, guessing a kind of normalisation?
e /= e.max()
plt.plot(e)
plt.show()

I know the whole len(e)/sr=time formula but I cant actually apply it to the graph. I try this:

plt.plot(len(y)/48000, e)

But I have np.shape issues:

ValueError: x and y must have same first dimension, but have shapes (1,) and (5846,)

How do I do this?

Edit:
I tried to do it using both of these (with error's afterwards):

plt.plot(y/48000, e)
ValueError: x and y must have same first dimension, but have shapes (2993006,) and (5846,)

plt.plot(e/48000, e)

Edit2:

plt.plot(np.arange(len(y))/48000, e)
ValueError: x and y must have same first dimension, but have shapes (2993006,) and (5846,)

plt.plot(np.arange(len(e))/48000, e) ended up with a graph with one vertical line down the middle.

-dear everyone this is what happens when you try do this with little background knowledge of the subject. Do your math homework!

You're getting an error because the two inputs to plt.plot need to be the same length. Notice that len(y)/48000 is just a number, so python tells us it of length 1.

Try plt.plot(range(len(e))/ 48000, e) **. What this does is change the values of your x axis to correspond to real time. If the sampling rate is $48000 \textrm{Hz}$, this means that $\textrm{sample}=1$ on your current axis now corresponds to 1/48000 of a second in the new axis, as this is the amount of time required to get one sample. Extend this to the whole axis - the time value on the new $x-$axis is simply the amount of time needed to take $y$ samples.

** if range() doesnt work, try numpy.arange() instead

Also, yes, these two lines

e -= e.min()
e /= e.max()

correspond to first shifting the signal so that its minimum value is now at 0, then dividing by the max. This is just a normalization, as you guessed.


Why do we use the top left corner as the origin in image processing? Why do we do this in image processing?:

Instead of this?:

I propose that the raster scan is a Western invention, in analogy with the way we write: left to right, top to bottom. Remember that "rastrum" denotes a  rake. When you rake a surface, or paint a floor, one usually does this such that one does not step on the raked or painted surface. The "raker" or "painter" moves backwards. One wants to avoid the situation described in Painted Into A Corner:

This defines a sort of 2D causality, using data for above and if not, from the same level (row) but left.  Parsing the 2D space top to bottom then left to right could have been an option. However, we know that the human vision is more sensitive to vertical than to horizontal, as one can see from the JPEG quantization matrix.

So perhaps, by going left to right first, the human mind could more easily spot a vertical change in a whole picture. *

However, the ancient slowness of a beam parsing a screen is not less important. For web browsing, it can be more important to first display a global low resolution patch and refine it iteratively, like in the progressive JPEG, JPEG 2000, or the clever Adam7 algorithm used in PNG.

See for instance http://www.libpng.org/pub/png/pngpics.html


audio analog input + usb data cable = sound noises I connect my android smartphone from audio jack output to my pc audio input and get a good analog sound signal. But when I connect this smartphone to  the pc with usb data cable the analog sound becomes hissing e.g. music from smartphone overlaps with noises on pc.
Why does this happen?
How can I eliminate these noises? maybe I have to use some shielded cables? If yes - then audio or usb? or I need some sound direct box or something like that?

The problem was resolved by adding ground loop noise filter like this


How do delay lines create so many effects? I was trialing this audio plug-in from SonicCharge (there's a program demo in the site).

It claims to be based on a delay line and uses some logic to modify it:

  At its core is a 12-bit digital delay with variable sample rate from $0$ to $352\textrm{ kHz}$. The delay is controlled by a programmable processor that allows you to change and modulate the delay time with various "operators".

I know from previous readings that delay lines can be used to construct a pitch shifter.

What are the properties in the delay line that allow for such numerous applications? I.e. what is it based on that the delay line can create such a varied amount of effects?

A DSP / Audio delay line doesnt create any effect in itself, a delay is simply a timed delay in the signal. it has no sound difference, it repeats the sound once and mixes it over original or sends it to a different mix line.

A feedback delay repeats the sound multiple times at the delay interval frequency, it sounds like a digital echo, and can become an oscillator above 20hz. almost all fx delays use variable feedback controls to control the number of repetitions of a sound, and they are filtered out as they feed back therefore diminishing.

Modern delays can work like grain synth samplers and often are based on them which have integrated effects like pitch shifting and research Grain-Delay and Grain-Resampling for infor about that. Squarepusher is somneone who uses alot of grain delay sometimes...

So a delay in effect is barealy an effect, but it can become an osciallator, has complex variants, and can feed into other lines that have any kind of additional FX.

To get a firmer grasp of that program i foudn this vid was more informative:

https://www.youtube.com/watch?v=ySzjFU6J4Ck

it's a pretty simple and very well balanced grain-delay with good controls and some bit crushing types of FX, and probably some filters too.


Extracting image metadata programatically I hope this is the right location to post this, as my question isn't really stack-overflow material but involves images.
I have hundreds of SEM-Images like this:

and need to track changes in certain features over the position at which the image was captured. This is why I would like to automatically retrieve metadata using something like bash or python for all images and save ti in a csv file or something.

The manufacturer of the SEM implements metadata on the images in the file itself, but I haven't had much success retrieving it. It doesn't show in the EXIF header if I inspect it with gimp or with image-magick's identify image.tiff. I only saw it if I save the image as an .xml file and view it in a texteditor.

The manufacturer ships an .exe program to inspect images, but as I am using Linux and have hundreds of images from which to retrieve the sample position, this isn't a very useful method for me.

Before I start writing a python / bash program that will convert the files to xml and use complex regex stuff to retrieve what I am looking for, I wanted to ask if anyone could imagine a simpler approach. My knowledge on image metadata is very limited.

Thanks in advance :)

Open a few files in a hex editor and check, if the header containing the information is always the same size, i.e. if the image data always begins at the same adress. If this is so, you can just retrieve the header with python ("normal" open, not binary). Then, you have the complete header as a string and can create any data structure you want, using any xml package you want.


K-means for 2D point clustering in python I have a set of points(2D) whose pixels are set and want to perform k-means on these pixels. Is clustering the 2D coordinates the right way ?

If so, can that be done using any libraries in python ?

It can be done very easily with the scikit-learn.  Examples are easy to find on their website, i.e. here. In my opinion it is the best way to go.

Modified code example from the above link:

import numpy as np
import matplotlib.pyplot as plt

from sklearn.cluster import KMeans
from sklearn.datasets.samples_generator import make_blobs

##############################################################################
# Generate sample data
np.random.seed(0)

batch_size = 45
centers = [[1, 1], [-1, -1], [1, -1]]
n_clusters = len(centers)
X, labels_true = make_blobs(n_samples=3000, centers=centers, cluster_std=0.7)

##############################################################################
# Compute clustering with Means

k_means = KMeans(init='k-means++', n_clusters=3, n_init=10)
k_means.fit(X)
k_means_labels = k_means.labels_
k_means_cluster_centers = k_means.cluster_centers_
k_means_labels_unique = np.unique(k_means_labels)

##############################################################################
# Plot result

colors = ['#4EACC5', '#FF9C34', '#4E9A06']
plt.figure()
plt.hold(True)
for k, col in zip(range(n_clusters), colors):
    my_members = k_means_labels == k
    cluster_center = k_means_cluster_centers[k]
    plt.plot(X[my_members, 0], X[my_members, 1], 'w',
            markerfacecolor=col, marker='.')
    plt.plot(cluster_center[0], cluster_center[1], 'o', markerfacecolor=col,
            markeredgecolor='k', markersize=6)
plt.title('KMeans')
plt.grid(True)
plt.show()

Yielding:


How to implement anisotropic diffusion using java? I am working on an image processing project in java where we use anisotropic diffusion for scale space generation.

In the below equation, after finding the diffusion value.using the second equation, I do not know how to find the divergence of diffusivity and gradient.

Since I am getting the diffusivity as a float value. I do not know how to find the divergence of that.can anyone please help me.

I coded a modest Java implementation of the Perona & Malik anisotropic diffusion, based on the ImgLib2 library. You can find it here.

You will have to import and learn how to use the ImgLib2 library though. Here is an example using the anisotropic diffusion algorithm.

Otherwise, here is a 2D only implementation, by Vladimir Pilny and Jiri Janacek, made for the ImageJ software.

And here is another one for ImageJ as well, by Christopher Mei.

Best regards.


OFDM Modulation: what's the output of the FFT? I'm trying to understand how OFDM is used in technologies like Wi-Fi.
On the Wikipedia page, an idealized transmitter is presented with this schema:

Data is converted from serial to parallel, mapped to some complex symbols using constellation mapping, then, using the inverse fast fourier transform, is converted to what seems one complex number with a real and imaginary part, which is then converted to some continuous signal using DAC.
On some other diagrams, the output of the IFFT is the same size as the input, and all of the outputs are multiplied by some (orthogonal) carriers before being transmitted, e.g. on this diagram (in French, but it should be clear what is what):

I'm having trouble understanding how these two diagrams relate with each other: are they different systems ?
As the FFT takes as an input a sampled continuous signal (so, let's say $N$ complex numbers), the IFFT should output $N$ complex numbers, but this does not appear on the Wikipedia diagram. Am I missing something?

not missing anything. You're right, if you have an N-point (I)FFT, the input is an N-element vector, and the output is an N-element vector.
The wikipedia diagram just "hides" the fact that these N output elements are "emitted" one after the other from the IFFT block.


Fourier transformations of discrete time signals How does one Fourier transform the following signals? $$x[n]=[0,9,0]$$ and $$y[n]=[9,0,19,0,9]$$
I've tried to get it along the way of $$X[\omega]=e^{-i\omega n}$$
but this seems incorrect. I'm not sure how the second one is supposed to be handled as well.

When finding the DTFT of short discrete sequences like this, you can directly compute the sum.

For example, say you had x[n] = [1 0 2] for n=0,1,2. Notice that I defined the values for n. If you had different n's, then the answer would be different.

You can compute the sum as
\begin{align}
X(\omega) &= \sum_n x[n] e^{-j\omega n} \\
& = 1e^{-j\omega 0} + 0e^{-j\omega 1} + 2e^{-j\omega 2} \\
& = 1 + 2e^{-j\omega 2}
\end{align}

Another note, you can sometimes simplify the answer by combining two exponentials. Be on the lookout for exponentials with the same coefficient. An example would be (using $e^{-j0}$ = 1)

\begin{align}
X(\omega) &= 2 + 2e^{-j\omega 4} \\
&= 2e^{-j\omega 0} + 2e^{-j\omega 4} \\
&= 2e^{-j\omega 2} (e^{j\omega2 } + e^{-j\omega 2}) \\
&= 2e^{-j\omega 2} (2\cos(\omega2)) \\
&= 4e^{-j\omega 2} \cos(\omega2)
\end{align}
Hope this helps.


How can I infer the cost function from Kruppa's simplified equations The following equations are Kruppa's simplified equations used in camera autocalibration. My objective here is to infer the cost function(Error Function) from this equations, So I can minimize the cost function through Levenberg-Marquardt nonlinear optimization algorithm.

(I already did the part of optimization, and only need the cost function in mathematical form so I can implement it).

I have read this paper, and this paper about this specific topic, but I didn't understand the complicated math inside.

Can anybody explain how the mathematical cost function will look like in simple language?

Thanks in Advance

$ W^{-1} = K K^{T}$,

where $W^{-1}$ is the dual of $IAC$, $K$ is the camera intrinsics matrix

$r,s$ are the first, and the second singular values of $S$

$u_1, u_2, u_3$ are the column-vectors of $U$,

$v_1, v_2, v_3$ are the column-vectors of $V$,

$USV = SVD(Fundemental\quad Matrix\__{ij})$,

and $i,j$ are two different images from moving camera.

According to the equation(30) in this paper The solution was to reform the Equations from the picture in my question here, and after that(Which the paper mentioned in different, complicated mathematical format) is to take the Frobenius norm to the subtract of newly formulated equations(30) on the paper.

$\rho_{1}/ \phi_{1} = \rho_{2}/ \phi_{2} = \rho_{2}/ \phi_{2}$ where $\rho$ is the numerator of each term of the equations, and $\phi$ is the denominator of each term, now the equation(30) in the paper proposed that: $\rho_{1}* \phi_{2} - \phi_{1}*  \rho_{2} = 0 $ and $\rho_{1}* \phi_{3} - \phi_{1}*  \rho_{3} = 0 $ now taking the Frobenius norm $\parallel.\parallel_{fro} $ for the two equations yields two values, Thus the cost function is the sum of this values.


How do I remove all white pixels in a modified image which are not present in the original? I have the following image

and after thresholding it, I obtained this result:

How can I remove all white pixels which do not belong to the original image?

The desired output should look like this (I manually edited the image):

Welcome to Image processing!

(See edit at the bottom)

If I understood you well, thresholding the image with a low threshold will provide a mask for the hand that you can use to remove the background pixels.

That's how it will look in Mathematica (ImagePad here is just used to remove the extra white border on your images):

img = ColorConvert[Import["http://i.stack.imgur.com/krbMs.jpg"], "Grayscale"];
border = BorderDimensions[img];
img = ImagePad[img, -border]

hand = Binarize[img, 0.2]

thresholded = Binarize[ImagePad[Import["http://i.stack.imgur.com/ylo3B.jpg"], -border], 0.1]

You may want to check this other question for an alternative segmentation technique based on ridge filtering.

EDIT: Based on your comment, you would like to only work on the thresholded image and clean it up to leave only the components corresponding to veins. I would look at morphological operations (e.g. Opening, components selection based on shape properties). By the way I feel this is harder than working on a different vein segmentation algorithm (see one proposed above).


time-domain channel estimation based on two vectors optimization Let's the input data vector $X = [X_1, X_2, X_3, X_4, X_5,X_6,X_7,X_8];$ where $[X_7,X_8]$ are well known, and the vector $y = h*X$ where $*$ is the convolution operation and $h = [h_1,h_2]$ is the channel vector.
I have a priory information that  $[X_7,X_8]$ = $T[X_1,X_2]$ where $T$ is a $2$x$2$ square matrix.
The issue I am trying to solve is can we get the channel vector $h$ based on optimizing $([y_7,y_8]$ - $T[y_1,y_2])$ ?  Is there any optimization algorithm we can use to estimate the vector $h$?

You specify that $T$ is known and invertible, so you know $X_1, X_2$, and then it's really trivial: $y[3]=[X_1,X_2]*[h_1, h_2]$, and $y[9]=[X_7, X_8] *[h_1,h_2]$; write that down as matrix system
$$\begin{pmatrix}y[3]\\y[9]\end{pmatrix}=\begin{pmatrix}X_1&X_2\\X_7&X_8\end{pmatrix}\begin{pmatrix}h_2\\h_1\end{pmatrix}$$
and solve it.


Question on convolution I am stuck at question number 2.8

This is how I have gone about solving it:

I have calculated $y(t)$ by convolving $x(t)$ with $h(t)$ using the fact that $x(t)$ convolved with an impulse at $t=t_0$ is just the same signal $x(t)$ time shifted by $t_0$. Hence,

$$y(t)= x(t+2) + 2x(t+1)$$

Now, according to the answer in the book, signal $y(t)$ is $t+3$ on $-2 < t \leq -1$. But, according to my answer, if we add $x(t+2) + 2x(t+1)$, the value of $y(t)$ at $t=-1$ is $4$. Why Please help me.

Also, why is it that in the answer in the book they haven't included $t=-2$ in the range?

$$x(t)=\begin{cases}
 t+1\qquad & 0\leq t\leq1 \\
 2-t\qquad &  1<t \leq2 \\
0\qquad & \mathrm{elsewhere}
\end{cases}$$

$$f_1(t)=x(t+2)=\begin{cases}
 t+2+1\qquad & 0\leq t+2\leq1 \\
 2-(t+2)\qquad &  1<t+2 \leq2 \\
0\qquad & \mathrm{elsewhere}
\end{cases}$$

$$f_1(t)=x(t+2)=\begin{cases}
 t+3\qquad &-2\leq t\leq-1 \\
 -t\qquad &  -1<t \leq0 \\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

$$f_2(t)=2x(t+1)=\begin{cases}
 2(t+1+1)\qquad & 0\leq t+1\leq1 \\
 2(2-(t+1))\qquad &  1<t+1 \leq2 \\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

$$f_2(t)=2x(t+1)=\begin{cases}
 2t+4\qquad & -1\leq t\leq0\\
 2-2t\qquad &  0<t \leq1 \\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

Now add them
$$y(t)=x(t+2)+2x(t+1)=\begin{cases}
 t+3\qquad & -2\leq t<-1 \\
3t+7\qquad &t=-1\\
 t+4\qquad &  -1<t \leq0 \\
2-2t\qquad &0<t\leq1\\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

Now check $f_1(-1)=2$,$f_2(-1)=2$,also $y(-1)=f_1(-1)+f_2(-1)=4$


Why is N-point DFT approximated by the sinc function? While looking into DFT leakage, I've came across the author saying that

  "..., the amplitude response of an N-point DFT bin in terms of the bin index m is approximated by the sinc function."

Other than the fact that the fourier transform of the rectangle function is the sinc function, I have no idea why the statement above is true.

When a waveform of arbitrary length is cut into a length that can be fed to an N-point DFT (with N being finite), that is the same as multiplying that original waveform by a rectangular window of that given finite length.  Multiplication in the time domain is the same as convolution in the frequency domain with the Fourier transform of that multiplying function, which happens to be a Sinc function for a rectangular window (or more precisely, for finite length real data, a conjugate mirrored Dirichlet (or periodic Sinc) function, which is pretty close to a regular Sinc, except near DC or Fs/2 or with very short DFTs.  More details here ).

Another way to look at it is that any pure sinusoid of a frequency that is not exactly equal to one of the DFT basis vectors will thus not be orthogonal to any of the DFT basis vectors.  Thus any non-zero magnitude of that sinusoid must end up partially represented in every element of a DFT result (since only orthogonal vector dot products will end up zero).  To interpolate the magnitude of that original non-orthogonal pure sinusoid, one might use the Whittaker–Shannon interpolation formula which happens to use a Sinc kernel to weight the contributions from all the DFT result elements into which energy from the original sinusoid was spread.


Aliasing below $f_s/2$ phi = exp(linspace(0, log(511), 1024)) - 1
x   = cos(2 * pi * phi)

Above will alias, despite peak instantaneous frequency evaluating to exp(log(511)) = 511 < 512. In fact, aliasing begins at x[724], long before peaking. The cause is revealed from diff(phi):

Successive inputs to cos differ by more than pi for all phi > .5. Thus while the continuous-time $\phi'(t) \leq 511$, the discrete-time greatly exceeds it, and we can't simply generate an exponential chirp via f0 * (k**t - 1) / log(k).
Can this limitation be predicted mathematically for any arbitrary $\phi(t)$ (more generally as in $\cos(\phi(t)$)? The goal is to know, for a given N, the allowed pairs of tmin and tmax in linspace(tmin, tmax, N) such that $\phi(t)$ won't alias. (-- Code)

My example isn't in fact below $f_s/2$; I blundered with $f_s=N$. Adjusting this eliminates aliasing, for linear, exponential, and hyperbolic chirps - and, in fact, for any signal with $\phi '(t) < f_s/2$.
This can be verified two ways:

For $\omega_\text{max} = \text{max}(\phi'(t))$ (assume $\phi(t)$ includes $2\pi$), the worst-case scenario is a straight line with slope $\omega_\text{max}$ (a curve can't go any steeper else it's no longer "max", and flatter is "more favorable"). But we know if $f_\text{max} < f_s/2$, a pure tone won't alias, thus the worst case is accounted for.

Aliasing occurs if the finite difference $\phi[n] - \phi[n-1] = \phi(n\Delta t) - \phi((n-1)\Delta t) \geq \pi$ for any $n$, where $\Delta t=$ sampling period, or increment of the t vector.  The greatest such difference is $\omega_\text{max} \cdot \Delta t$ - setting the inequality: $\omega_\text{max} \Delta t < \pi \Rightarrow \omega_\text{max} < \pi / \Delta t \Rightarrow f_\text{max} < 1/(2\Delta t)= f_s/2$.

Thus if I ensured $f_s/2$ was met, phi could never exceed $\pi$.

Example in question, corrected: if we seek (tmin, tmax) = (0, log(511)), our $N$ must be large enough so that $\Delta t$ is small enough for $\text{max}(\phi'(t)) < \pi$. The max is 2*pi*exp(log(511)) = 2*pi*511 = wmax, so max finite difference is wmax*dt < pi -> dt < pi/wmax, and dt = t[1] - t[0] = (tmax - tmin)/N, thus (tmax - tmin)/N < pi/wmax -> N > wmax/pi*(tmax - tmin) = 2*pi*511 / pi * log(511) = 2*511*log(511) = 6373.6.
N = 6374
t = np.linspace(0, np.log(511), N, endpoint=False)
phi = 2*np.pi * (np.exp(t) - 1)
x   = np.cos(phi)


Differences between matlab's 'impulse' and 'lsim' command; create impulse response I've got following problem:

assume you've this continous system 1st order:

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

s = tf('s');
sys = 1/(1+s);

I'm going to discretize this system via impulse invariant transformation and plot the results.

sysD = c2d(sys,1e-3,'impulse');
impulse(sys);
hold all;
impulse(sysD);

So this seems to work.

If I try to create the impulse response for 7 seconds manually I use:

imp = zeros(1,7*1e3);
imp(1) = 1;
t = 0:1e-3:7-1e-3;
y = lsim(sys,imp,t);
plot(t,y)

Unfortunately this is not the response I expected. The gain seems to be scaled by the samplerate 1e3.

So the question is: What does impulse do since this response should be correct. lsim uses the filter command and simply passes the nominator and denominator to filter I think. Is the dirac impulse (or kronecker delta) wrong? Must it be scaled by the samplerate?

I discovered the whole thing, because I wanted to design a very simple digital filter. I knew the responses of time continuous systems, so I thought it would be easy to discretize the transfer function and use the coefficients of the digital transfer function that can be used by filter or in a dsp application.

What is wrong?

Thank you very much!

Indeed the lsim function is returning scaled values. You might also notice that it is connected to the gain $k$ of your filter:

[z,p,k] = tf2zpk(sysD.num{1}, sysD.den{1});

That will return:

k=0.002

You can scale it back by the sampling period as you noticed.


Kalman filter for tracking sinusoidal motion I am wanting to create a Kalman filter that can be used to track an object undergoing sinusoidal (lets assume simple harmonic) motion. I have seen many examples and implemented my own python code for a simple SUVAT equation/projectile motion case however I have been having trouble finding examples where a Kalman filter is used to track an object moving with sinusoidal motion and I have been having trouble working out how to construct the State transition matrix (commonly called A) and the Control matrix (commonly called B) to model such behaviour.

This will eventually be used to track a noisy voltage signal corresponding to the position of an oscillating object such that one can recover it's 'true' motion.

I would appreciate if someone could direct me to some useful resources where I could learn how to do this or explain how I might go about doing this.

Assume that the amplitude remains constant as well as the angular frequency $\omega$. The phase will be predicted using $\phi_k=\phi_{k-1} + \omega T_s$

The measurement matrix usually denoted as H and is  the Jacobian matrix or the first derivatives matrix of the model with respect to the state variables. B will be the all zeros matrix.


Averaging out noise and SNR improvement Suppose I have a sequence of identical symbols carrying the same information at a sample rate equal to the symbol rate, why do we always talk about averaging them to average out noise, why not just  adding them all together without  averaging them? I suppose the improvement in SNR is going to be the same assuming that noise is uncorrelated from one sample or symbol to another, is this related to the limited dynamic range in DSP?

Yes, they both give you the same SNR gain. I suppose there might be a dynamic range issue, as you say, if a LOT of samples are added, but I think the main reason that the literature speaks of "averaging" instead of "adding" is because for systems that allow a variable integration rate, we can be assured that other parts of the system work the same, like AGC.


Books that explain DSP well to those not directly in engineering? I do work with computer graphics and am dipping my toes into ray tracing.  That field involves a good number of the subjects covered in DSP (Fourier transform, time vs frequency space, etc) but I was hoping to find a book in the realm that didn't also deal with, say, digital devices or electrical engineering-- not that I have anything against those subjects, but they're not germane to what I'm doing and I'd prefer not to wade through them on my way to DSP understanding.

Could anyone recommend some books for someone in my position?  I've taken calculus and differential equations & linear algebra in the past and so am reasonably proficient as far as the mathematics involved.  Thanks much.

I am surprised no one has mentioned Richard Lyon's book - by far one of the BEST books out there on understanding digital signal processing in a very clear, concise, and methodological way.

Its excellence comes from that fact that he explains concepts to you in a very easy way to grasp, without loss of rigor or detail needed to get to the heart of various DSP concepts. He does not shy away from math, (which is good), but you are guaranteed to understand all of it.


Control systems and convolution I think i am not understanding the concept of convolution well.

Lets say we are given a system impulse response in the S-domain, and we have implemented a controller $G_c(s)$ that will adjust the system response to some kind of input. Now when we are implementing this controller say using an operational amp. circuit all we need to do is to get the differential equation of $G_c(s)$ and start implementing it with the correct component values. Is this correct?

And if so, lets say we are trying to implement this controller on a computer, where our inputs, outputs are no longer continuous signals. And our controller in the Z-domain is represented by $G_c(z)$. Now the next step to implement this controller in a computer is to get the difference equation and implement it in a software. Is this correct?

And if so, what is the difference between implementing the controller difference equation on a software and implementing the convolution algorithm on the software?

For example if we are given g[n] as an array of samples representing the the controller unit impulse response, and our input is stored inside x[n] a simple for loop might do this algorithm for us. As far as i know

for i in len(x):      #loop through the input elements
 for j in len(g):     #loop through the controller unit impulse response elements
  y[i+j] += x[i]*g[j]  #Store the conv. result in the output array y

Which is just telling the software to scale every point in the impulse response by x[n] and then shift the response by n, then add them all up.

I am probably missing something here !

This's a too broad question including continuous time systems, discrete time systems, control systems, convolution operation, difference operation and differentiaiton operations.

However your main focus is on convolution. Therefore I will state in plain simple terms that the operation of convolution defines and calculates the output of a linear time invariant system to a given input from the impulse response of the system. This is performed by the convolution integral for continuous time systems and convolution sum for discrete time systems. Therefore convolution forms the core of signal processing. Just as integration or differentiation forms the core of calculus.

Control theory is an application of extensive signal processing, feedback theory, linear algebra, calculus, matrix theory, transform theory and probability theory. Therefore it will naturaly utilize convolutions in its many formulations. Nevertheless convolution is nothing more than a tool for the control field point of view.

A physical process requires an analog (continuous time) system to control it. But modern technology enables digital computers (sampled data systems) to be used instead. As a consequence, integrations and differentiations of the continuous time system are replaced by sums and differences in digital software implementation of analog systems. Whereas the link between the digital and continuous time systems can be defined via sampling theorem or some other mapping relations.


Shifting a fourier spectrum by subpixel amount in python I am working on a Fourier Ptychography problem. My research problem requires me to shift the Fourier spectrum of an image by a floating-point value. For a real-valued image, we can simply use cv2.warpAffine to shift the image by floating-point values (aka subpixel shifting).

Taking fourier transform of an image produces a complex matrix. The problem is, cv2.warpAffine does not support complex matrices, and so I cannot use it on them. I tried searching for alternatives, but none of them seem to work. I came across numpy.roll, but the problem is, it does not support subpixel shifting. Rounding off the shifting values translates to loss of information in my case. Is there a solution in python, that allows for subpixel shifting on complex matrices?

Thanks.

EDIT: Based on Marcus' answer, I did some digging and implemented a nifty little script for subpixel shifting in python based on a Matlab script for the same.

Here's the link to the script. Hope it helps!

You'd not do it in frequency domain at all. Use the shift property of the discrete Fourier transform!

Simply multiply the rows of your original image with an $e^{j2\pi \Delta f_x x/W}$ pointwise ($x$: pixel index in that row, $W$: width of image) before transforming to frequency domain to shift by $\Delta f_x$ in row direction. Same for shifts in column direction; multiply columns with $e^{j2\pi \Delta f_y y/H}$ pointwise ($y$: pixel index in that column, $H$: height of image). There's no restrictions on the "fineness" of $\Delta f$ in either direction, and any 2D shift can be understood as a shift in row and one in column direction.

If necessary, transform your image from frequency to spatial domain, do the multiplications, and transform back.


Proving conditions for controllability Let's say I have the following LTI system:
$$\dot{x}(t) = \mathbf{A} x(t) + \mathbf{B} u(t)$$
I need to somehow show the following is true or false (proof):

This system is controllable if and only if for any $\mu$ such that $(−\mu\mathbf{I − A})$ is Hurwitz, there exists a unique positive deﬁnite solution $\mathbf W$ to the Lyapunov equation:

$$\mathbf A \mathbf W + \mathbf W \mathbf A^T − \mathbf B\mathbf B^T = −2\mu\mathbf W$$

Can anyone see the proof? I tried replacing $\mathbf A$ in the Lyapunov equation with $(−\mu\mathbf I − \mathbf A)$, and then adding the $−2\mu\mathbf W$, to the left side, but the $\mu$'s end up cancelling out, so I'm not sure.
I'd really appreciate any help.

Not really a direct help, but this paper might provide some insight, particularly the paragraph below.


What is the difference between Antenna port and Physical Antenna in LTE? I have just started working on the 5G NR physical layer and I have difficulty understanding the concept of Antenna Port at the downlink.

In release 10 there supposedly 15-22 antenna ports for CSIRS (Channel State Information Reference Signal) and different transmission modes have different antenna ports assigned to them how are these antenna ports actually mapped to physical antennas.

As I understand each port has a resource grid for itself how does this antenna ports resource grid map onto a physical antenna?

Is it possible to have two antenna ports but only one physical antenna or one antenna port but two physical antenna's? If yes then how does the Resource elements map onto the physical antenna's?

This answer is for LTE standard, not 5G NR. However I do believe that the basic concept does not change since.

(Logical) Antenna port is defined formally as in 3GPP 36.211

  An antenna port is defined such that the channel over which a symbol on the antenna port is conveyed can be inferred
  from the channel over which another symbol on the same antenna port is conveyed.

This simply means that symbols that are transmitted over an antenna port are
subject to the same propagation conditions. Hence to determine the characteristic channel for an antenna port, a separate channel estimation for each antenna port needs to be done. That is the reason why an antenna port is characterized by its Reference Signal.

To the best of my knowledge, each antenna port has its own resource grid that is mapped directly to the corresponding physical antenna port. Then if two antenna ports are mapped to the same physical antenna port, their resource grid are summed up. There will be reserved resource elements to avoid inteference between antenna ports.

The four first antenna ports (1-4) are normally mapped 1-1 to physical anntenna ports. And it is possible to have two antenna ports mapped to one physical antenna port. For example, antenna port 5 can be mapped to the same physical antenna ports of (logical) antenna ports 1-4.

It is also possible that one (logical) antenna port is mapped to two (or more) physical antenna ports. It is the case of beamforming that requires several physical antennas to create a beam that is characterized by UE-RS. The beamforming effect is transparent to UE and UE only sees one resource grid estimated by its UE-RS.

This white paper can provide more useful details.


Matlab IIR cheby2 bandpass, problems with Magnitude I am trying to implement an IIR bandpass in Matlab. There are two things I don't understand. But first of all, let me post the code.

fs = 50000; %Hz
f1filter = 1000;  %Hz
f2filter = 2000;  %Hz
filterOrdnung = 6;   %erzeugt (2 * filterOrdnung) Filterkoeffizienten
anzahlFilter = 5;    %Anzahl der zu berechnenden Filter zur Auswertung
schrittweiteVerschiebung = 1000;  %Verschiebung in Hz
astop = 20;           %Verstärkung an den Bandgrenzen des Bandpasses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k = zeros(anzahlFilter, 1);                 %Vektor für Gain

z = zeros(2*filterOrdnung,anzahlFilter);    %Matrix für
                                            %Nullstellen

p = zeros(2*filterOrdnung, anzahlFilter);   %Matrix für
                                            %Polstellen

sos = zeros(anzahlFilter * filterOrdnung + (anzahlFilter - 1), 6);  %Matrix für Filterkoeffizienten
[z(:,1), p(:,1), k(1,1)] = cheby2(filterOrdnung, astop, [((2*f1filter)/fs)...
    ((2*f2filter)/fs)]);

[sos(1:6,:), g] = zp2sos(z(:,1),p(:,1),k(1,1));    %Output form:[b01,b11,b21,1,a11,a21
                                                   %             b02,b12,b22,1,a12,a22
                                                   %             b0N,b1N...........a2N]

fvtool(sos(1:6, :),'Analysis','freq');

So I want to create a bandpass of the bandwith of 1000Hz to 2000Hz. The filter order is the 6th order and the stop attenuation should be 20dB. I want to have the Structure of second order form.
(By the way, is there any specific Matlab function to create an IIR filter as a cascade of Directform II structures?)

I am creating "sos" much larger than I need it to be now because I want to create several different bandpasses, but let's just focus on this example.

1)First question: Why do I receive 12 numerator and 12 denominator coefficients when I only have 6th order IIR? In the books I have read, that for a 6th order IIR I should receive 7 numerator and 7 denominator coefficients, with the first denominator coefficient = 1.

2)At the header of the code, I declare astop as 20dB. The Matlab documentation of cheby2 tells me, that this is the attenuation of the stopband. But when I let the code run and look at the magnitude, the magnitude of the passband is about 20dB. I expected it to be 0dB while the stopband should be -20dB. The normed frequencys seem to be about correct though.
I hope I could explain my problems well enough.

Actually we are also facing the same problem, instead of (Z,P,k), even if you use (A,B,C,D), as soon as sos will come in syntax with g, your magnitude response will scaled by the attenuation factor. If you are concerned only with plotting magnitude response you can plot as freqz(sos) you will get exact magnitude response without scaling.


Cascade Transfer Function; Odd or Even Transfer Function suppose we have a cascade realization of a transfer function of a higher order.

What difference does it make if the order of the transfer function is an Even or an Odd number?

If the order $N$ of the system is even, you get $N/2$ second-order sections, where each second-order section combines either a pair of complex conjugate poles and zeros, or two real-valued poles and zeros. If $N$ is odd, you get $(N-1)/2$ second-order sections and one additional first-order section for the remaining real-valued pole and zero.


Using MATLAB to define a function in signal processing 101 I want to define a function $y(t)$ that has 2 cases:

$$y(t)=
\begin{cases}
-1 & \text{for} & -T \lt t \leq 0 \\
1 & \text{for} & 0 \lt t\leq T
\end{cases}
$$
where $T$ is an integer.

I already defined $t$ using linspace between $-T$ to $T$ in a resolution of 1000 points.

I was advised to use the command ones, but I dont know how to apply it by the function terms.

Thank you.

Well, this is Octave but should work in Matlab:

function y=Y(t)
  y=ones(size(t));
  y=y+(t<=0)*(-2);

.. and it is vectorised. Calling will look like:

>t=linspace(-T,T,1000);
>y=Y(t);

To make this more general we could also make the time at which the function switches a parameter also:

function y=Y(t,S)
  y=ones(size(t));
  y=y+(t<=S)*(-2);


Technique to detect pitch after analysing only one pitch period? does anyone know if there is a technic to estimate the pitch of a given signal by analysing only one pitch period.

I'm asking because I read about a program called AudioMotors.
https://lesound.io/product/audiomotors-pro/
If I'm right, you can load in a recording of your car, where you recorded the rpm from idle constantly rising to max rpm. Then you can jump to any position u want and hold the rpm at this point. There are a lot more features this program provides. But I'm only interested in the part where you can hold a chosen rpm and let it sound pretty good.

I'm not an expert in dsp but I thought you can jump to a position you want. Then you select a few samples. From these samples you have to estimate the pitch of the first period and match the pitch of all following pre selected samples (pitch scaling). After this you can loop these samples and not hear any pitch change. And maybe let it sound ok.

If you can help me with this or have any suggestions on how to program a "light" version of AudioMotors pls let me know.

Thanks in advance.

in order to estimate the period of a quasiperiodic signal, normally you need an entire period plus a little more.  you need to compare a little snippet of the waveform with its counterpart in the adjacent period.  that's the fastest i expect a pitch detector to come up with anything.


Usefulness of Matched $z$ transform Method I'm aware that the matched $z$ transform method maps between the continuous $s$ plane and the discrete/digital $z$ plane but my question is - when would this be necessary? Why would we need to convert between the two?
Thanks

the Matched Z method is the simplest method to convert an analog filter design with transfer function $H_\mathrm{a}(s)$ to a digital filter design $H(z)$.  It does it by mapping every pole and every zero from the $s$-plane to the $z$-plane using:

$$ z \leftarrow e^{sT} $$

where $T \triangleq \frac{1}{f_\mathrm{s}}$ is the sampling period and the reciprocal of the sample rate, $f_\mathrm{s}$.  Stable analog filters get mapped to stable digital filters.  That's all that it is.

The other two common methods are the Bilinear Transform, which approximates the above mapping:

$$ z \leftarrow \frac{1+\frac{sT}{2}}{1-\frac{sT}{2}} \approx \frac{e^{\frac{sT}{2}}}{e^{\frac{-sT}{2}}} = e^{sT} $$

and the Impulse Invariant method which is a time-domain way of looking at things:

$$ h[n] \leftarrow h_\mathrm{a}(nT) $$


Discrete time inverse fourier transform of cosine squared $$ X(\omega) = \cos^2(\omega)$$

I tried this problem, and I ended up getting $0$, which doesn't make any sense. I integrated:

$$ x(n) = \frac{1}{2\pi}\int_{0}^{2\pi} \cos^2(\omega)e^{j{\pi}n} d\omega $$

Which can be changed to

$$ x(n) = \frac{1}{2\pi}\int_{0}^{2\pi} \bigg(\frac{e^{j\omega}+e^{-j\omega}}{2}\bigg)^2e^{j\omega} d\omega $$

To
$$ x(n) = \frac{1}{2\pi}\int_{0}^{2\pi} \frac{e^{j\omega(n+2)} + e^{j\omega(n-2)} + 2e^{jn\omega}}{4} d\omega $$

Now, taking out the four, each anti-derivative term contains the original term. And since each original term is of the form $$e^{jk\omega}, k\in \mathbb{Z}$$ and we're integrating from $0$ to $2\pi$, which are identical frequencies, the result ends up being zero. What went wrong?

The only thing that went wrong is your last statement:

They're only zero if $k\ne 0$.

Thus, you get three Dirac delta distributions.

Let's look at your last equation, and work with that:

$$
\begin{align}
 x(n) &= \frac{1}{2\pi}\int\limits_{0}^{2\pi} \frac{e^{j\omega(n+2)} + e^{j\omega(n-2)} + 2e^{jn\omega}}{4}\,d\omega \tag1\label{orig}\\
&=\frac1{8\pi}\left(\int\limits_{0}^{2\pi}e^{j(n+2)\omega}\,d\omega + \int\limits_{0}^{2\pi}e^{j(n-2)\omega}\,d\omega + \int\limits_{0}^{2\pi}e^{jn\omega}\,d\omega \right)\tag2\label{linearint}
\end{align}$$

From $\eqref{orig}$ to $\eqref{linearint}$ I just went ahead and used the fact that integration is a linear operation to decompose that single integral into a sum of three, very similar, integrals.

Now, let's look at the last one, as you (sadly, with a mistake), already did:
$$
i(n):=\int\limits_0^{2\pi}e^{jn\omega}\,d\omega =
\begin{cases}\frac1{jn}e^{jn\omega}\Big|_{\omega=0}^{2\pi} \overset{e^x\text{ is }j2\pi\text{-periodic}}= 0& \text{for } n\omega \ne 0\\
\int\limits_0^{2\pi}e^{0}\,d\omega= \int\limits_0^{2\pi}1\,d\omega= 2\pi& \text{else.}
\end{cases}
$$

So, this integral yields zero for all values of $n$ but for $n=0$, in which case it yields a fixed, finite, non-zero value.

Your other two integrals are just the same, but shifted w.r.t. $n$ by $2$ to the left or the right.

For more info on the Dirac delta distribution (looks like a function, isn't quite exactly a function), I'd recommend your favourite textbook or wikipedia. It's a thing that, when integrated, has an integral of $1$, but in an arbitrarily narrow support around $0$, outside of which its value is $0$.


Very simple question about signal periodicity $$x[n] = u[n]+u[-n]$$

Is it periodic or not? My answer is

$$u[n] = {1       , n\geqslant0}$$
$$u[-n] = {1      , n\leqslant0}$$

which means that the signal $x[n]$ is always equal to $1$ from $-\infty$ to $+\infty$, but it equals $2$ at $n=0$. Am I right?

no, but $x[n]$ is, what we call, an "even-symmetry" function:

$$ x[-n] = x[n] \quad \forall n \in \mathbb{Z} \ . $$

what makes $x[n]$ periodic is $$ x[n+N] = x[n] \quad \forall n \in \mathbb{Z} $$ which is really the premise for the DFT of length $N$.


Fisher ratio or F_ratio I am confused while writing code to implement Fisher’s F-ratio.

Suppose I have 4 speakers, each speaker pronounces 2 words and each word has 3 frames and each frame has 30 sub-bands.

$u$ is mean of all speakers.

Suppose $x_{ij}$ is one sub-band energy of the $j$th speech frame of
speaker $i$ with $j= 1\ldots N$, and $i = 1 \ldots M$. $u_i$ and
$u$ are the sub-band energy averages for speaker $i$ and for
all speakers, respectively, which are defined as
$$
u = \frac{1}{MN} \sum_{i=1}^{M}\sum_{j=1}^N x_{ij}\\
u_i = \frac{1}{N} \sum_{j=1}^N x_{ij}
$$
I have written following MATLAB code:

U=mean(mean(mean(mean(spectr1(:,:,:,:)))));
for subband=1:N
nominator=0;
for speaker=1:speaker_number
   nominator=(mean(mean(spectr1(speaker,:,:,subband)))-U).^2+nominator;
end
denaminator=0;
denaminator=(mean(mean(mean(spectr1(:,:,:,subband))))-mean(mean(spectr1(speaker,:,:,subband)))).^2;
F_ratio(subband)=((1/speaker_number)*(nominator))/((1/(speaker_number*phonem_number))*denaminator);
end

Well, check out Wikipedia.

The ratio is defined as:

$$
F = \frac{\mathrm{between-group\ variability}}{\mathrm{within-group\ variability}}
$$

where, using your notation, they say that the between-group variability is:

$$
\sum_{i=1}^{M} (u_i - u)^2/(M -1)
$$

and the within-group variability is:

$$
\sum_{i=1}^{M} \sum_{j=1}^{N}( x_{ij} - u_i )^2 / (N - M)
$$


Appropriate sampling frequency for bpsk modulated signal I'm modulating a message with particular baudrate, say 2MHz using BPSK , Carrier frequency(Fc) is 8MHz, What should be appropriate Sampling Frequency(Fs) following Nyquist criteria, Should it be greater than 2*(2*baudrate) or 2*(Fc+baudrate)?

The sampling frequency depends not only on the carrier frequency and the symbol (baud) rate, but also on the pulse shape. If you transmit with square pulses, you'll need a very large bandwidth (in theory, infinite). If you transmit with sinc pulses, then the base-band bandwidth is half the pulse rate. Other pulses will fall somewhere in between.

To expand on Peter K.'s comment: you should look into using the complex envelope for simulation. The advantage is that, when you use it, the sampling frequency depends only on the symbol rate and pulse shape, not on the carrier frequency.


What is the tightest known bound on the reconstruction error in compressed sensing? I am specifically thinking about the reconstruction error of L1-minimization from compressed measurements with noise.
I know a bound from (8) in The restricted isometry property and its implications for compressed sensing, but I was wondering if there is anything newer and or tighter that I am not aware of. Additionally, I would like to relate $\delta_{2s}$ to the size of e.g. a Gaussian measurement matrix. It seems I can do that from Lemma 3.1 and (3.22) in Decoding by linear programming, but how tight is (3.22) and can I get closer?

According to Michael McCoy on LinkedIn:

The best bound for Gaussian matrices is in: W. Xu and B. Hassibi, Precise Stability Phase Transitions for l1 Minimization: A Unified Geometric Framework, IEEE Trans. Inf. Theory. Oct. 2011

See, for example, the curves in Fig. 1.


O(1) peak segmentation algorithm in matlab? I have an intensity image, similar to the image below.  I would like to very quickly, segment this image, into separate peaks.  Where each peak is labelled, and segmented.  A notional segmentation of 3 blobs is given below by the black lines.

Definition of a peak segment is:  The contour that just crosses the highest saddle point between this peak and an adjacent peak.

This has to be a well studied problem, but I can't figure it out.  I have been looking at the Matlab functions:

imregionalmax
imreconstruct
imhmax

all of which are close.... It is not hard to write an algorithm to do this (see below) using imreconstruct to operate on the tallest peak, segment out the first peak, and repeat.  The algorithm below is O(N) in number of peaks.  Such an algorithm can and should be O(1) in number of peaks (and this is important for my application.)

The above image was rendered in matlab with the command

imagesc(peaks)

Here is an algorithm that performs the segmentation:

a = peaks;
peakValsX = [21, 25, 35];
peakValsY = [20, 38, 25];

N = numel(peakValsX);

clear('peakImgs', 'peakMask');
for ix = 1 : N;
    aptImg       = zeros(size(a));
    aptImg(peakValsY(ix), peakValsX(ix)) = a(peakValsY(ix),peakValsX(ix));
    peakImgs(:,:,ix) =  imreconstruct(aptImg,a);
    % All we've done is eliminated the other peaks.
    assert(all(reshape(peakImgs(:,:,ix) <= a,[],1)));
end

for ix = 1 : N
    peakMask(:,:,ix) = ix*(max(peakImgs(:,:,setdiff(1:N, ix)),[],3) < a);
end

peakMaskResult = max(peakMask,[],3);
imagesc(peakMaskResult);

And the resulting segmentation:

To Reiterate.  The question is, "Is there a matlab implementation of the peak segmentation algorithm (as described above) that can be executed in O(1) of number of peaks?"

I'm not sure why you need to call imreconstruct more than once. The way you normally use morphological reconstruction is described here: You basically subtract some constant (the peak threshold) from your source image, then reconstruct it. What you get is an image that is identical to your original image except at the peaks:

So you simply take the difference and binarize, or compare the two images pixel-wise.


Values of angular frequency of discrete-time sinusoidal signal Why are the values of angular frequency of discrete-time sinusoidal signals in between -pi and pi?

Because the sequence

$$ x[n] = A \cos( \omega n) $$

will produce exactly the same samples for all discrete-time frequencies such as $\omega = \pi/3$ or $\omega = \pi/3 + 2\pi$ or $\omega = \pi/3 + 6\pi$ or $\omega = \pi/3 -8 \pi$... In other words every frequency of the form $w = w_0 + 2\pi k $ for all integers $k$ will yield the same set of samples; they are indistinguishable, and we chose to use the range $[0,2\pi]$ or $[-\pi,\pi]$, when referring to the discrete-time sinusoidal frequencies.

Note that the periodicity of the discrete-time frequency $\omega$ relies on the fact that the index $n$ is an integer. Otrherwise, like in the continuous-time case, the frequency $\omega$ would range from $-\infty$ to $\infty$.


Convex Optimization in Signal and Image Processing In signal processing, convex optimization plays a useful role in problems such as sparse signal recovery and filter design. What other places does convex optimization appear?

For example, in compressed sensing the Basis pursuit Denoising problem, the LASSO problem and the Dantzig selector can be posed as:

\begin{eqnarray}
\min_{x} \ell(Ax-b)+r(x)
\end{eqnarray}

where $\ell(\cdot)$ is and $r(\cdot)$ are appropriate loss and regularization terms, respectively. Moreover, the design of a filter subject to time and frequency constraints often yields a convex formulation.

There's a whole area of signal processing dedicated to optimal filtering. In pretty much every case I've seen the filtering problem is formulated with a convex cost function.

Here's a freely available book on the subject - Sophocles J. Orfanidis - Optimum
Signal Processing.


Inverse Fourier Transform of a squared sinc function I have here a squared sinc function, which is the Fourier Transform of some triangular pulse:
$$\mathrm H(f)= 2\mathrm A\mathrm T_\mathrm o \frac{\sin^2(2\pi f \mathrm T_\mathrm o)}{(2\pi f\mathrm T_\mathrm o)^2}$$
As an excercise, I would like to go back to the original time domain triangular pulse, using the inverse Fourier Transform.
I feel like I'm very close to achieving it, however, I stumbled upon this integral:
$$\mathrm h(t)=\frac {\mathrm A} {2\pi ^2 \mathrm T_\mathrm o}\int_{0}^{\infty}\frac{(1-\cos(4\pi f\mathrm T_\mathrm o))\cos(2\pi ft)}{f^2}df$$
I'm not sure if there is already an expression which defines the value of this integral, or if this integral is of a well-known type, with a well-defined value. The tables of integrals that I have right now are basic and don't cover integrals of such type.
Does someone know if there is a well-known expression for the value of integrals of this type?
Or, would you be so kind as to recommend a good table of integrals for Fourier Analysis?
Thank you very much!

For $T_o>0$, this integral can be tackled with a Schwinger parameterization and Fubini's Theorem.
(The equivalent of the intermediate integral you specifically asked about is at the 7th equal sign below, with the trig functions written in terms of complex exponentials.)
$$\begin{align*}h(t) &= 2AT_o\int_{-\infty}^{\infty} \frac{\sin^2(2\pi fT_o)}{(2\pi fT_o)^2}e^{2\pi ift} \; df\\
\\
&= 2AT_o\int_{-\infty}^{\infty} \frac{\sin^2(2\pi fT_o)}{(2\pi fT_o)^2}\cos(2\pi ft) \; df\\
\\
&= 4AT_o\int_{0
}^{\infty} \frac{\sin^2(2\pi fT_o)}{(2\pi fT_o)^2}\cos(2\pi ft) \; df\\
\\
&= 4AT_o\int_{0}^{\infty} \sin^2(2\pi fT_o)\cos(2\pi  ft)\int_0^{\infty}ue^{-u2\pi fT_o} \; du\; df \\
\\
&= 4AT_o\int_{0}^{\infty} \int_{0}^{\infty} u\sin^2(2\pi fT_o)\cos(2\pi  ft)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{i2\pi fT_o}-e^{-i2\pi fT_o}\right)^2\left(e^{i2\pi  ft}+e^{-i2\pi  ft}\right)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{i2\pi f 2T_o}-2+e^{-i2\pi f 2T_o}\right)\left(e^{i2\pi  ft}+e^{-i2\pi  ft}\right)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{i2\pi f (2T_o+t)}+e^{i2\pi f (2T_o-t)}-2\left[e^{i2\pi  ft}+e^{-i2\pi  ft}\right]+e^{-i2\pi f (2T_o-t)}+e^{-i2\pi f (2T_o+t)}\right)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{-f(2\pi uT_o - i2\pi  [2T_o+t])}+e^{-f(2\pi uT_o - i2\pi  [2T_o-t])}-2\left[e^{-f(2\pi uT_o -i2\pi  t)}+e^{-f(2\pi u T_o+i2\pi  t)}\right]+e^{-f(2\pi uT_o + i2\pi  [2T_o+t])}+e^{-f(2\pi uT_o + i2\pi  [2T_o-t])}\right) \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty}  u\left(\dfrac{1}{2\pi uT_o - i2\pi  [2T_o+t]}+\dfrac{1}{2\pi uT_o - i2\pi  [2T_o-t]}-\dfrac{2}{2\pi uT_o -i2\pi  t}-\dfrac{2}{2\pi u T_o+i2\pi  t}+\dfrac{1}{2\pi uT_o + i2\pi  [2T_o+t]}+\dfrac{1}{2\pi uT_o + i2\pi  [2T_o-t]}\right) \; du \\
\\
&= -\dfrac{A}{4\pi}\int_{0}^{\infty}  \dfrac{u}{u - i\left[2+\frac{t}{T_o}\right]}+\dfrac{u}{u - i\left[2-\frac{t}{T_o}\right]}-\dfrac{2u}{u -i\frac{t}{T_o}}-\dfrac{2u}{u+i\frac{t}{T_o}}+\dfrac{u}{u+ i\left[2+\frac{t}{T_o}\right]}+\dfrac{u}{u + i\left[2-\frac{t}{T_o}\right]} \; du \\
&= -\dfrac{A}{2\pi}\int_{0}^{\infty}  \dfrac{u^2}{u^2 +\left[2+\frac{t}{T_o}\right]^2}+\dfrac{u^2}{u^2 + \left[2-\frac{t}{T_o}\right]^2}-\dfrac{2u^2}{u^2 +\left[\frac{t}{T_o}\right]^2}\; du \\
\\
&= -\dfrac{A}{2\pi}\int_{0}^{\infty}  \dfrac{u^2}{u^2 +\left[\frac{t}{T_o}+2\right]^2}+\dfrac{u^2}{u^2 + \left[\frac{t}{T_o}-2\right]^2}-\dfrac{2u^2}{u^2 +\left[\frac{t}{T_o}\right]^2}\; du \\
\\
&= -\dfrac{A}{2\pi} \left(u - \left[\frac{t}{T_o}+2\right] \tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}+2}\right]+u - \left[\frac{t}{T_o}-2\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}-2}\right]-2u +2\left[\frac{t}{T_o}\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}}\right]\right)\bigg{|}_0^\infty\\
\\
&= \dfrac{A}{2\pi} \left(\left[\frac{t}{T_o}+2\right] \tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}+2}\right]+\left[\frac{t}{T_o}-2\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}-2}\right]- 2\left[\frac{t}{T_o}\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}}\right]\right)\bigg{|}_0^\infty\\
\\
&= \begin{cases} \dfrac{A}{2\pi} \left(-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]+2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right) & \dfrac{t}{T_o} \le -2 \\
\\
\dfrac{A}{2\pi} \left(\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]+2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right) &-2 \le \dfrac{t}{T_o}\le 0  \\
\\
\dfrac{A}{2\pi} \left(\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]-2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right)&0 \le\dfrac{t}{T_o}\le 2 \\
\\
\dfrac{A}{2\pi} \left(\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]+\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]-2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right)&\dfrac{t}{T_o}\ge 2\\
\end{cases} \\
\\
h(t) &= \begin{cases} 0 & \dfrac{t}{T_o} \le -2 \\
\\
A\left(1+\dfrac{1}{2}\cdot\dfrac{t}{T_o}\right) &-2 \le \dfrac{t}{T_o}\le 0  \\
\\
A\left(1-\dfrac{1}{2}\cdot\dfrac{t}{T_o}\right)&0 \le\dfrac{t}{T_o}\le 2 \\
\\
0 &\dfrac{t}{T_o}\ge 2\\
\end{cases}\\
\\
\end{align*}$$