What is the impulse response used in an interpolation filter when upsampling? I need to downsample a single period waveform from $M$ samples to $L$ samples. If that matters, in my particular case, $L=M-1$.

I read about downsampling and came to the conclusion I need:

expand the original waveform by a factor $L$ (some other source call the step “zero stuffing”),
interpolate the missing samples,
decimate by a factor of $M$:

If I understood it correctly, expansion only requires adding $L-1$ zeros between each original sample. Decimation seems simple too since it only requires keeping every $M^{th}$ samples.

But I'm stuck with the interpolation process. The naive solution of using linear interpolation gave me surprisingly good results at low frequencies but quickly degraded when the frequency increases (i.e., .when the number of samples in the waveform decrease).

I found several mentions of using an interpolation filter that looks like a dot product1. For example, on the Upsampling Wikipedia page, I see:

  $$y[j+nL] = \sum_{k=0}^K x[n-k]\cdot h[j+kL],\ \ j = 0,1,\ldots,L-1,$$

  where the h[•] sequence is the impulse response, and K is the largest value of k for which h[j + kL] is non-zero

In my understanding, $x[•]$ is the sequence after the expansion step. $y[•]$ is the upsampled waveform. But I don't understand what are $h[•]$ and $K$.  I read they are related to the impulse response and I understand the concept. But I can't see where it comes from and how to obtain it in practice.

Could someone explain to me how to obtain $h[•]$ and $K$?

1It also looks like a convolution to me, but I'm not quite sure of that. And apparently this would be an instance of FIR filter. We already talked about that in other questions, but this is still not clear to me.

The dot in that summation is just scalar multiplication.  And yes, it's a convolution -- you're convolving the input signal by the filter.


A cubic interpolation function: folkloric copypasta or clever trade-off? I've been reading on interpolation methods recently and I have come across an implementation of cubic interpolation that is leaving my head scratching. Every other variant and example of cubic interpolation I've come across has included lots of fractions and exponents. This implementation I found on Paul Bourke's website uses only additions and subtractions to figure the coefficients. I've seen this implementation used in several codebases but unfortunately without any documentation.
double Cubic(double x, double a, double b, double c, double d) {
    double A = d - c - a + b;
    double B = a - b - A;
    double C = c - a;
    double D = b;
    return A * (x * x * x) +
           B * (x * x) +
           C * x +
           D;
}

If I may, let me first walk through how this implementation might have come to be just to make sure I at least understand the algebra going on here. But I'm curious as to why it can be made so simple.
As I understand it, cubic interpolation takes four sample points, say $a$, $b$, $c$, and $d$, and interpolates between b and c using the function:
$$
f(x) = Ax^3 + Bx^2 + Cx + D \tag{1}\label{1}
$$
Where $0 \le x \le 1$.
A major condition I often see (one that's obviously important for me) is that $f(x)$ must run through the sample points. So, in this case, samples $b$ and $c$ correspond to $x = 0$ and $x = 1$ respectively.
$$
f(0) = b \tag{2}\label{2}
$$
$$
f(1) = c \tag{3}\label{3}
$$
There is another condition that says that the derivative at samples $b$ and $c$ are the same as the slope between their surrounding points. I believe this is called a Catmull-Rom spline?
$$
f'(x) = 3Ax^2 + 2Bx + C \tag{4}\label{4}
$$
$$
f'(0) = \frac{c - a}{2} \tag{5}\label{5}
$$
$$
f'(1) = \frac{d - b}{2} \tag{6}\label{6}
$$
Now, if I follow that path, I get Paul Breeuwsma's solution here and that all makes sense to me.
However, in order to get Bourke's implementation, I have to multiply the slopes at $b$ and $c$ by 2! That is:
$$
\begin{align}
f'(0) = c - a \tag{5a}\label{5b}\\
f'(1) = d - b \tag{6a}\label{6b}
\end{align}
$$
If I follow from that, I get the below for the coefficients $A$, $B$, $C$, and $D$ in order to match Bourke's.
First, $D$:
$$
\require{cancel}
$$
$$
\begin{align}
f(0) &= \cancel{A(0)^3} + \cancel{B(0)^2} + \cancel{C(0)} + D \\
f(0) &= D \\
D = b \tag{7}\label{7}
\end{align}
$$
Now, $C$:
$$
\begin{align}
f'(0) &= \cancel{3A(0)^2} + \cancel{2B(0)} + C \\
f'(0) &= C \\
C = c - a \tag{8}\label{8}
\end{align}
$$
Now, $B$:
$$
\begin{align}
f(1) &= A(1)^3 + B(1)^2 + C(1) + D \\
f(1) &= A + B + C + D \\
c &= A + B + c - a + b \\
B = a - b - A \tag{9}\label{9}
\end{align}
$$
Finally, $A$:
$$
\begin{align}
f'(1) &= 3A(1)^2 + 2B(1) + C \\
d - b &= 3A + 2(a - b - A) + c - a \\
d - b &= 3A + 2a - 2b - 2A + c - a \\
d - b &= A + a - 2b + c \\
A = d - c + -a + b \tag{10}\label{10}
\end{align}
$$
Ok, so... why the heck? This looks like nothing I've seen while reading up on this.
By doubling the slope at $b$ and $c$, besides having a fast and definitely appealing implementation, is this a trade-off or are there additional conditions that were able to be made that I missed? It seems like this could cause the interpolant to shoot around much more but I haven't yet coded up a comparison to look at and listen to. I've been scouring articles, blog posts, and papers and I can't seem to map anything to this. If anything, other implementations end up quite hairy and difficult for me to understand. Is there a name for this kind of interpolation? Where did it come from?
I know this was a long question. Thanks for reading!

A reasonable, analytically derived cubic interpolation method will give a straight line if the inputs are on a straight line. With $a = 0$, $b = 1$, $c = 2$, $d = 3$ you get from cubic_bourke $A = d - c - a + b = 3 - 2 - 0 + 1 = 2$ as the coefficient for $x^3$ whereas for a straight line you would have $A = 0$ (and $B = 0$ as the coefficient for $x^2$). So the method doesn't give a straight line for straight-line input. Instead it will, as you say, wobble. In the frequency response that will manifest itself as a low multiplicity of the roots at multiples of the sampling frequency.
The impulse response is not symmetrical. Normally splines have a symmetrical impulse response. The asymmetry will result in phase non-linearity:
Figure 1. Impulse response of cubic_bourke.
The magnitude frequency response shows a rather nice treble response for audio purposes, but gives a lot of aliasing, worse than linear interpolation:
Figure 2. Magnitude frequency rseponse of cubic_bourke (red), linear interpolation (blue) and cubic Hermite spline (green). Half the input sampling frequency is at $\omega = \pi$.
Cubic Hermite (Catmull-Rom) will give better treble response up to 15 kHz for a 44.1 kHz input sampling frequency, and lower aliasing, than cubic_bourke.
So I say scrap it.


"Bi Directional" Kalman Filter I am working on a project in Object Tracking, i.e. need to predict the location of next bounding box.

I used a Hungarian algorithm with a Kalman Filter (which is a common method in this domain) which produced decent results. However, lots of the times detector has False Positive or False Negative detections (i.e. noise in the sensor) which limits Kalman's predictions.

I was thinking to apply Kalman filter in the backward time direction, similar to bidirectional recurrent neural networks. This way it would utilize information in the future, of course compromising on some latency. This way we would have 2 filters, one in forward direction and another in backward both of which would predict a current location and two results will be merged using another simple method. From my understanding of Kalman filters it would only be possible to have a backward filter only if I were to recompute the full state at every sequence step given next n states.

Anuar Y, Welcome to the DSP community.

What you're talking about is called smoothing.
Let me explain, assume we have samples $ {\left\{ x \left[ n \right] \right\}}_{n = 0}^{N - 1} $ and we want to build estimator for $ x \left[ k \right] $ which we will define as $ \hat{x} \left[ k \right] $.

Now, we have 3 types of estimation:

The case $ k > N - 1 $ is called prediction.
This is usually what we employ Kalman Filter for.
The case $ k = N - 1 $ is called filtration.
This is easy as we basically filter data we have.
The case $ k < N - 1 $ is called smoothing.
This is basically what you're after.

So, there is plenty of information about the Kalman Filter in the Smoothing framework:

Derivation of Kalman Filtering and Smoothing Equations.
Kalman Filtering and Smoothing.
Fixed Leg Smoother.

If there is something specific you need, let us know and we'll assist you.


Does the Shannon theorem not apply when the amplitude of a wave is changed faster than half the time period of the wave? Shannon's version of the sampling theorem states that if a function contains frequencies all strictly less than $B$ hertz, then it is completely determined by giving its ordinates at a series of points spaced $\frac{1}{2B}$ seconds apart.
Now, let us suppose we are talking about just one frequency (say a laser with frequency $B-\varepsilon$ hertz, where $\epsilon$ is an arbitrarily small positive number) with an amplitude that is varied over time for relaying a signal. Let's assume that samples of this wave are taken at the following points of time, measured in seconds:
$ \{0, \frac{1}{2B}, \frac{2}{2B}, \frac{3}{2B}, \ldots \} $.
Suppose the laser amplitude is changed to half its value at $\frac{1.4}{2B}$ seconds and is then doubled back to its original amplitude at $\frac{1.6}{2B}$, wouldn't this information go undetected? It seems to me that the theorem assumes that amplitude does not change inside one wavelet.

Once you start changing the amplitude you are increasing the bandwidth of the signal. That's called "amplitude modulation" and the highest frequency is now the sum of the original frequency and the highest frequency in the modulation signal.
The sampling theorem still holds. You still only need twice the bandwidth but the bandwidth has increased due to the modulation


GnuRadio CC Decoder or How to use FEC Extended Decoder I am currently trying to recover a satellite signal, encoded following CCSDS standards. This includes a convolutional code with rate 1/2 and constraint length 7.
I am receiving the signal through an SDR directly connected to the sender using coax cables and some attenuation in-between. For testing purposes, I am sending 0 and 1 alternating. For demodulation and decoding, I am currently using the graph attached to this post.
Up until the convolutional decoder, everything seems to work without problems. When deactivating the code both on the sender and receiver, I can decode the data without any problems. When the convolutional code is active, the decoder sometimes is able to decode the data properly, but most of the time, the output just contains garbage.
I therefore wonder, am I using the decoder block in a wrong way? Documentation about how to properly use this block is hard to find. I already tried different settings for the decoder definition but so far nothing helped.
Can anyone help me with this?
Kind regards
valkyrie

Inside CCSDS channel coding standard, there is a synchronization word called the ASM.
It is a pattern sequence of 32, 64 or 256 bits which is used for synchronization inside the stream. It is usefull to solve the two points I describe below.
Inside your GNU block diagram I see two things that might be an issue:

Phase ambiguity: There is no mechanism to solve the phase ambiguity of the BPSK in demodulation. I am not really sure if Viterbi decoder works fine if not solved (the ambiguity).

Two bits ambiguity of CC: This type of code with a rate of 1/2 produces two encoded bits for one noncoded. At the input of decoder, you must be synchronized with the first encoded bit of a couple.

This could explain that sometimes your system works after switching on/off. On chance over two or maybe four...
Other remark, using an alternating pattern of 0 and 1 could make you miss something debuging. Try a 32 bits pattern.
Tell us if you found something


Given a Horizontal Derivative Kernel What Would Be the Output Image? Hi, im solving this problem and everything went well until i went and solved the position p(5,3) and i got the  answer of 0. However, the solution guide  lists the solution  as 150 which is pointed at in blue. I don't understand this since the pixel at the location is 0, meaning we don't have to calculate the solution however the solution says otherwise. Thanks for your time.

The solution provided is correct. Here is why.

To compute the derivative at index(3,5)

d(3,5) = value_at(4,5) - value_at(2,5)
d(3,5) = 150 - 0
d(3,5) = 150

Hence the solution provided there is correct. Hope this helps :)


Metric Spaces: Why $L_\infty$ selects the maximum value I have a basic question about the metric spaces. There are several metric spaces like $L_1$, $L_2$ to $L_\infty$. The $L_p$ metric is defined by the following equation:

$$d_p(x,y)=\left(\sum_{i=1}^{n}|x_i-y_i|^p \right)^{1/p}$$

Now, when $p = \infty$, then the equation of the matrix is following:

$$d_{\infty}(x,y)=\max_{i=1,2,\ldots, n}|x_i-y_i|$$

Now, I am taking the value of $x$ as $x = [1,2,3]^T$ (Transpose) and compute $L_p$ metric for, $p = 1,2$ and $\infty$.

My question is, why $L_\infty$ metric choose the maximum value.

In the general case, let  $x = (x_1,\dots,x_n)$ be a finite-length vector (in a finite dimensional space).
The finite sequence of absolute values $|x_i|$ does attain its maximum (because the sequence is finite), denoted $M = \max_i |x_i|$.

Let $m$ be the (exact) number of coordinates in $x = (x_1,\dots,x_n)$ whose absolute value  is equal to $M$. Thus, $1\le m\le n$. Then, we can lower and upper bound the $\ell_p$ norm of $x$ as follows: $$(mM^p)^{\frac{1}{p}}\le\ell_p(x)\le (nM^p)^{\frac{1}{p}}\,.$$

Both $m^\frac{1}{p}$ and $n^\frac{1}{p}$ tend to $1$ as $p\to\infty$, thus $\ell_p(x) \to M$, by the squeeze theorem.

[EDIT] To provide more concrete substance, let us see what happen with your example: $x=[1,2,3]^T$ (the transposition does not change the result): $$\|x\|_p = d_p(x,0) = (1^p+2^p+3^p)^{1/p} = 3\times\left(\left(\frac{1}{3}\right)^p+\left(\frac{2}{3}\right)^p+1\right)^{1/p}\,.$$

Both $\left(\frac{1}{3}\right)^p$ and $\left(\frac{2}{3}\right)^p$ tend to $0$ as $p \to \infty$. Thus, $\|x\|_p \to 3$.


What is the total impulse response in a system with feedback interconnection? Let's assume that we have a system with a typical feedback interconnection where the output is given by the following equation:

$$y(t) = \left(x(t) - z(t)\right) \star h_{1}(t) \tag{1} $$

where:
$$z(t) = y(t) \star h_{2}(t) \tag{2}$$

So, is $h_{1}(t)$ the system's total impulse response or the convolution between $h_{1}(t)$ and $h_{2}(t)$ after replacing $z(t)$ in the equation $(1)$ while using the associative property?

It is easier to work in the $s$-domain:

$$Z=YH_2$$
$$Y=(X-Z)H_1$$
Hence,
$$Y=(X-YH_2)H_1=XH_1-YH_1H_2\Rightarrow Y(1+H_1H_2)=XH_1$$
Therefore,
$$H(s)=\frac{Y(s)}{X(s)}=\frac{H_1(s)}{1+H_1(s)H_2(s)}$$
which is called the closed-loop transfer function. The closed-loop impulse response can be found by inverse Laplace transform.


Baby Cry Sound detection I wanted to write a code to detect baby cry sound.
I am using Windows as platform.  Presently, I am able to get audio samples and its frequency plot(using FFT) but not sure how to proceed forward.

I wanted to ask what steps I should follow to detect the baby cry sound given its time-frequency plot.

I saw some methods such as median filter followed by HMM in speech recognition. But for simple sound detection do I need to go for such sophiticated method?

I will be very grateful if you could help me.

Since you have frequency bucket information, you essentially have a frequency histogram.

If you treat histograms as n dimensional vectors, you can use dot product to get a similarity value that goes from 0 to 1 and could tune that so that sounds that have a value greater than say 0.6 are considered a match.

For more information check this out:  http://blog.demofox.org/2015/02/13/writing-a-basic-search-engine-aka-calculating-similarity-of-histograms-with-dot-product/

Alternately you could get a little more advanced and use support vector machines, which require training data of both negatives and positives to find a way to differentiate hits and misses.


Linear Equalizers Weights For linear equalizers, the weights of the equalizer are set to achieve minimum error between the equalized signal and the original input signal. First, a training sequence is sent as input, and the weights that achieve the minimum error between the equalized and the original training sequence are obtained. My question is:

In case of non-adaptive equalizers, are the weights obtained for the training sequence used for any other input signals?
In case of adaptive equalizers, how the weights are updated according to the new input signals? The error should be minimized between the equalized signal and what? There is no known sequence in this case.

The error should be minimized between the equalized signal and what? There is no known sequence in this case.

But there is: if the channel is still good enough so you can decode a packet successfully, you can reconstruct the original signal by simply re-encoding the packet.


How to Smooth Gradient Estimates for Steepest Descent Optimization In steepest descent
methods of minimizing a function $f(x), x \in \mathbb{R}^d$,

it's common to approximate the gradient by finite differences:
$\qquad\qquad \nabla f(x) \approx gradest( x; h ) \equiv
{ {f( x + h I ) \ - \ f(x)} \over h } $
then minimize along those directions:
$\qquad\qquad x_{k+1} \  \text{from min} \ f( x_k + \lambda \ gradest( x_k ) ) $

Can a signal-processing point of view suggest ways of smoothing $gradest()$ in this context ?
For example, one could view the "zig-zags" in the picture above as high frequencies in $gradest()$,
to be smoothed out over the last few ${x_k}$.
But that may be naive:
can one model, then smooth, zig-zags through non-uniformaly spaced points in 5d or 10d ?

In the Probabilistic settings we have many methods applied to the Stochastic Gradient Descent in order to decrease the variance of the Gradient Estimation (ADAM / RMS Prop / AdaDelta, etc...).

The nice thing is to utilize them in deterministic settings.
So for instance you can use Momentum which to Signal Processing guy will look just like applying IIR / AR filter on the Gradient.

There are also methods derived from the Convex Optimization world settings.
There is nothing prevents us using them in any setting besides there are no guarantees.
For instance you can use FISTA Method or Nesterov Accelerated Gradient Descent.

More References:

Deniz Yuret's Homepage  - Alec Radford's animations for optimization algorithms.
Sebastian Ruder - An overview of gradient descent optimization algorithms.
Shao-An Lu - SGD > Adam?? Which One Is The Best Optimizer: Dogs-VS-Cats Toy Experiment.
FAST AI (Sylvain Gugger and Jeremy Howard) - AdamW and Super-convergence is now the fastest way to train neural nets.


efficient spectrum/fft analysis of sparse signals (found in radar applications) consider the classic scenario of identifying a radar frequency. You typically have a very short pulse of very high RF frequency. For example, a 50ns pulse repeating itself every 1MHz or so (PRF 1MHz and Pulse Repitition Interval 1us) and let's say there are 50 such pulses in each (radar) interrogation so that each set of coherent pulses takes 50us.

The RF frequency is generally in Gigahertz but let's assume it has been downconverted to lie somewhere between 200MHz to 400MHz. Now you want to analyze this signal after subsampling it at a frequency of 100MHz. The straightforward approach would be to sample the whole 50 us worth of data but in reality the only time is signal is on during the 50us period is 50 times 50 ns which is just 2.5 us. The rest of the signal is just zeros. My question therefore is:

is there a 'smart' way of getting rid of this apparently useless information and hence be able to do a much shorter fft?
is it possible to throw away long strings of zeros in between active pulses, stick the pulses together and perform fft? would that cause any sort of distortion in the spectrum?

One solution for this kind of problem is called an Instantaneous Frequency Measurement (IFM) Receiver. The pulses are located with an envelope detector and a threshold. The average phase change per sample is computed across each pulse; this is proportional to frequency of the pulse.


Autocorrelation and PSD Let $X(t)$ and $Y(t)$ be two orthogonal processes with power spectral densities
$$S_{xx}(f) = S_{yy}(f)=\begin{cases}
1-\lvert f\rvert, & \lvert f\rvert<1
\\[1ex]
0,& \text{otherwise}
\end{cases}$$

Define a new process $Z(t) = Y(t)-X(t-1)$. Determine and sketch the power spectral density $S_{zz}(f)$.

The solution is the following:
The autocorrelation function of the process is
$$R_{zz}= E \bigg\{\big[Y(t-\tau)-X(t+\tau-1)\big]\big[Y(t)-X(t-1)\big] \bigg\} $$
And now comes the part which I do not understand at all.
\begin{align}
R_{zz}&=R_{yy}(\tau)-R_{yx}(\tau+1)-R_{xy}(\tau-1)+R_{xx}(\tau)\\
&=R_{yy}(\tau) + R_{xx}(\tau)\ \quad\text{since}\quad R_{yx} = R_{xy} = 0 \quad\text{from orthogonality.}
\end{align}

Therefore, $S_{zz}(f)=2S_{yy}(f)=2S_{xx}(f)$ as shown below.

Can anyone explain it like I am five?

HINT:

Expand the terms inside the expectation $E\{\ldots\}$, and then use the fact that the expectation of a sum equals the sum of the individual expectations (due to linearity).


How to detect rectangular images inserted into another image? It's just for a simple binary classification, so it's not needed to find the location of the inserted image or the content. I'm only interested in checking whether something like that is present or not.
Positive examples:

Negative examples:

The inserted image rectangles are always upright, i.e., they are not rotated by non-90-degree steps. They cover at least 1 % of the image area, i.e., their edges are at least 10 % of the image dimensions.
The images have undergone lossy compression to an unknown degree.

I've tried looking for straight horizontal and vertical edges using Scharr filters, Canny edge detection, and Hough lines, but the results are not very reliable. Any creative ideas would be appreciated. :)

My current attempt looks as follows:
import numpy as np
from cv2 import cv2

def get_value(image: np.ndarray) -> float:
    assert len(image.shape) and image.shape[2] == 3, "please provide a color image"

    scharr_scale = 1 / 16

    # Detect vertical and horizontal edges on all three channels (BGR) separately.
    scharr_x_abs = np.abs(cv2.Scharr(image, cv2.CV_32F, dx=1, dy=0, scale=scharr_scale))
    scharr_y_abs = np.abs(cv2.Scharr(image, cv2.CV_32F, dx=0, dy=1, scale=scharr_scale))

    # Reduce edge images to one channel.
    scharr_x_abs_gray = cv2.cvtColor(scharr_x_abs, cv2.COLOR_BGR2GRAY)
    scharr_y_abs_gray = cv2.cvtColor(scharr_y_abs, cv2.COLOR_BGR2GRAY)

    # Erode the edge images perpendicular to the filter direction to remove noise.
    scharr_x_abs_eroded = cv2.erode(scharr_x_abs_gray, np.ones((5, 1)))
    scharr_y_abs_eroded = cv2.erode(scharr_y_abs_gray, np.ones((5, 1)))

    # Get the average edge value per column/row part.
    scharr_x_row_abs = cv2.resize(scharr_x_abs_eroded, (scharr_x_abs_eroded.shape[1], 8), interpolation=cv2.INTER_AREA)
    scharr_y_col_abs = cv2.resize(scharr_y_abs_eroded, (8, scharr_x_abs_eroded.shape[0]), interpolation=cv2.INTER_AREA)

    # Build second derivative, but perpendicular to the original one,
    # to detect abrupt starts and ends of edges.
    scharr_x_row_abs_scharred_again = np.abs(cv2.Scharr(scharr_x_row_abs, cv2.CV_32F, dx=1, dy=0, scale=scharr_scale))
    scharr_y_col_abs_scharred_again = np.abs(cv2.Scharr(scharr_y_col_abs, cv2.CV_32F, dx=0, dy=1, scale=scharr_scale))

    # The maximum abrupt start/end of edges are the final result.
    scharr_x_row_max = np.max(scharr_x_row_abs_scharred_again)
    scharr_y_col_max = np.max(scharr_y_col_abs_scharred_again)
    result = float(np.maximum(scharr_x_row_max, scharr_y_col_max))

    return result

Output for the provides examples images:

100.99776458740234: A high value for this positive example. That's good.
37.71356201171875: A low value for this positive example. That's not good.
56.305137634277344: A not-too-high value for this positive example. That's ok.
27.330982208251953: A not-too-high value for this positive example. That's ok.

So there is no threshold that would divide the two classes correctly.

In the two positive examples, the rectangular box is "likely" to be noticed visually, notably by the frame, but also the difference of the content inside and outside the box, especially in the color diversity.
However, the compression (if lossy, JPEG_like) may induce local fake edges (quantization of the DC and higher frequency), local noise and fake colors along $8\times 8$ squares.
Hence, I would suggest a local/global algorithms in three steps.

First, estimation: a crude segmentation with closed contours. For instance with geodesic active regions or level sets. The principle is as follows: from a starting contour, the border evolve until it meet a region that either has a gradient (edge) and/or the statistics of the regions on both sides are different enough. Generally they are performered on colored gradient-tensors. Although the contour has some regularity, it is allowed to split and separate into several region. I think the is doable by maximizing the color contrast and allowing sharp edge
Second: model matching. You get a set of connected contours, which can be fitted by Cartesian rectangular boxes. The latter can be probed with your size constraints. One question remains: how to select only one box?
Third: decision/machine learning. Given a set of images, and a training ensemble of negatives and positives examples, you can train a logistic classifier on features related to the inner/outer boxes.
If you cannot learn examples, a distance between the feature-vector of the inner box to the feature-vector of the outer box could provide a quantitative metric related to some probability of having a rectangular image  inserted.


Modal analysis and frequency response function - interpretation needed For couple of weeks now, I have been trying to measure cutting forces during machining processes. The system which I am currently examining is highly distorted at high frequencies due to its structural dynamics. My set-up consists of force sensor, workpiece holder and workpiece itself - they are screwed all together with M12 screw.

The problem is that I struggle to interpret FRF of the setup correctly. For better understanding, please see attached figure.

What is the physical explaination of the peaks marked in the figure? Does it mean that forces measured at these frequencies are damped?

Thanks in advance for any kind of help.

Best,
Daniel

This may be a better question for https://physics.stackexchange.com/ but I'll give it a shot.

I assume you show the FFT of the time waveform at the sensor. In order to get the actual frequency response you would need to divide by the FFT of the hammer response. Since the hammer is reasonably smooth, we can eyeball this. It also looks like your measurement is "good": sufficient signal to noise ratio and nothing unexpected in there.
Force is a vector: it's unclear whether you are only measuring the normal force or a full 3 component vector and calculate the magnitude.
At a broad brush, the frequency response can be interpreted as "reasonably flat" with a bunch of local wiggles at specific frequencies
The flat part of the frequency response (FR) simply means the piece itself is moving (or compressing) uniformly: i.e. it's moving as one without the sensor surface deforming.
The wiggles are modal behavior. Here the piece (or surface) has different movement at different locations.
The shape of the modes is determined by the geometry (and material) of the piece. That's why the modes tend to be very narrow band.
Not all modes will be visible in a measurement like this. Modes have "nodes" and "peaks" as a function of position. In order to see mode both excitation and sensor need to be far enough away from a node. If the hammer is close to the node, the mode will not be excited and if the sensor is close to the node, the sensor won't pick it up (despite it being there and active).
Modal movement implies different areas of "positive" and "negative" movements.
The modal movement simply adds to the "move as one piece" movement. If both go in the same direction you see a peak and if both go in opposite directions you see a dip.

  Does it mean that forces measured at these frequencies are damped?

Not really: what it means is that movement due to modal behavior and movement of the entire piece as a single chunk are in opposite directions at this particular point and frequency. So the sum is smaller than each individual movement.


Cannot recreate a simple Nyquist series in Matlab a Nyquist series should look something like {...-1, 1, -1, 1, ... }, which i think means it is a sine wave sampled at the peaks and valleys exactly.

Now, given a sampling rate of 48 Hz, the Nyquist frequency should be 24 Hz.

If I am not mistaken, Nyquist series is basically the the dataset of a sine wave of frequency = Nyquist Frequency, sampled at the sampling rate given above.

in short, these are my matlab vars:

SR = 48;             //sampling rate (hz)
X  = 0:1/SR:2*pi;    //domain (time in general)
fNy= SR/2;           //nyquist frequency
Y1 = sin(2*pi*fNy*X);// a sine wave with freq = nyquist frequency.

when i plot X and Y1, i get a really messed up plot. (i think its zero since its so small in magnitude)

I believe that Y1 should be exactly = 1, -1, 1, -1, 1, -1 ...
Where did I go wrong? thanks

HINT:

Try Y1 = cos(2*pi*fNy*X);

(And think about the $x$-axis values at which you sampled the sine)


Downsampling lowpass filter for audio: FIR or IIR? I'm working on a real-time audio processing project in which I have to downsample a 44.1 kHz to a yet-to-be-determined lower sampling rate essentially for reducing computational complexity. I'm looking for a decent downsampling filter, and I found articles suggesting both IIR (mainly elliptic) and FIR (mainly minimum phase). FIR minumum phase filter wouldn't necessarily work for me, because I need a linear phase filter. So the options are either linear phase FIR that could have a considerable latency in it or a double-filtering IIR (MATLAB filtfilt command) that first filters forward then backward.

What are some pros and cons of FIR and IIR as downsampling filters? Which is more practical in a real-time application?

If linear phase is a requirement, that will probably steer you toward an FIR implementation. It is possible to build IIR filters that have approximate linear phase, but it is easy to design a linear-phase FIR.

If you're concerned about latency, forward-backward filtering as in filtfilt isn't really a good option. In general, it's really meant to be used an offline process, since to implement the technique exactly, you have to run the entire signal through forward, then do the same in reverse. That implies that you have access to the entire signal at once, which is not commensurate with low delay.

In general, an FIR filter will require a higher order for a given set of performance requierments. However, FIR filters bring some real advantages, such as guaranteed stability, lower susceptibility to roundoff errors (since the quantization error doesn't get fed back through the filter, although you can compensate for this with some increased complexity), and simply-achieved linear phase response. In addition, efficient FIR filter implementations are available for many processor architectures, mitigating the cost of the extra taps somewhat.

Another way to mitigate the extra-tap cost for FIR filters in your situation is to take advantage of efficient multirate signal processing techniques. Specifically, you can use a polyphase decimation approach to significantly reduce the number of computations that you perform in the process of decimating your signal. This has the effect of decreasing the number of effective taps (in terms of computational complexity) in the decimation filter. In addition, if you need to decimate by a large factor, then multistage approaches can help further reduce your load. Lyons' introductory DSP book has some good easy-to-read material on these topics.

Given more specific parameters of your system, one could make more pointed recommendations. What are your filter design requirements? What sort of computational capabilities does your platform have? What sample rate will you decimate to?


Well known waveforms or signal equations for filter testing I'm testing a software based DSP filter. To test I'm generating some sine waves and adding them together. I'm about to start trying to generate more complicated input signals with trigonometric functions. My question is whether there is a set of well known waveforms with algorithms that can be programmatically generated?

A simple example I currently use to generate an input waveform is:

yl_sum =  np.sin(2*np.pi*x / (freq * 3)) + np.cos(4*np.pi*x/(freq*2)) + np.cos(7*np.pi*x/(freq*2))

I was thinking there might be an existing set of algorithms along these lines that are a good basis to start testing with?

My question is whether there is a set of well known waveforms with algorithms that can be programmatically generated?

Yes and no. There are many wave forms that can be created programatically. However these may or may not be useful for you test.

Your test should focus on determining whether the algorithm meets the stated requirements. Hence the best way to build a test is to run down the list of specifications and requirements and build test that clearly establish whether the implementation meets the requirement or not . Once you know what exactly a specific test needs to assess and what the pass/fair criteria are, you can start thinking about how to build the test. That includes what wave forms but also what test equipment, what test flow, what data acquisition, what metrics etc.

Wave forms are a means to an end. Once you have a clear definition of what the end is, the wave form question will most likely become obvious.


partial derivative of image I have to find the partial derivative of an image with respect to its x dimension. I am using central difference method i.e

$ \partial_x F(x) = \frac{F(x+1,y)-F(x-1,y)}{2} $

Here $F(x,y)$ represents the image and if I want to use spatial filtering for the same then I can use filter mask as

$0.5 \times [0, -1, 0;0, 0, 0;0, 1, 0],$

and then simply do filtering of F with this mask.

Now I want to know:

whether central difference can be applied for finding partial derivative?
whether the filter mask (which I have proposed above) can be used or not for finding partial derivative?
if not, please suggest me a filter mask which can perform it better and also mention how is it performing the required operation.

Computing the image derivatives are very important for a lot of vision tasks. So let me first begin with answering your questions quickly and then get to more details.

Whether central difference can be applied for finding partial derivative?
YES
Whether the filter mask (which I have proposed above) can be used or not for finding partial derivative?
It can be used, but only the vertical one (in $y$ direction). The horizontal one would be the transpose of that.
If not, please suggest me a filter mask which can perform it better and also mention how is it performing the required operation.

For noisy images computing derivatives this way would not be very helpful. That's why Sobel operator exists (http://en.wikipedia.org/wiki/Sobel_operator). Instead of the mask you proposed, it suggests to use $0.125 \times [1, 0, -1; 2, 0, -2; 1, 0, -1]$. However Sobel lacked the rotational symmetry. So Scharr kernels were proposed (the same wiki page, on the very bottom). These were proven to be better alternatives. Yet, these developments did not stop the further research. Being not so old, Farid and Simoncelli proposed more robust differentiation operators. See the paper:

Hany Farid and Eero Simoncelli "Differentiation of Discrete
Multi-Dimensional Signals" IEEE Trans. Image Processing. 13(4): 496-508 (2004)

One implementation of this is presented in Peter's Computer Vision MATLAB functions:

http://www.csse.uwa.edu.au/~pk/research/matlabfns/

Check derivative5 and derivative7 functions.

Reading about Scale Space would also give you more insight on the nature of such linear operators: http://en.wikipedia.org/wiki/Scale_space

Hope these help. Cheers,


Audio: How to measure dynamic range (in sense of dynamic compression/ "loudness war") I'm about to calculate some quality measures for audio.
One important of them is the dynamic range (in sense of "loudness war").
By which in this context is not meant the SNR but
the (naive said) "variance" of the "occuring energy".
https://en.wikipedia.org/wiki/Loudness_war
http://dr.loudness-war.info/
My idea was to take the average RMS for that and use
its standard deviation (relative to the average RMS).
But this doesn't separate the good from the bad.
If I take TT DR Offline meter it outputs DR11 for
a fine (Steven Wilson) remaster,
and DR7 for a deadly compressed one from an everyday compilation.
But the relative standard deviation of them is nearly
the same (about 69% and 65% for L and R, window = 50ms).
Are there any suggestions/ ideas how to measure it better?

One important of them is the dynamic range (in sense of "loudness war"). By which in this context is not meant the SNR but the (naive said) "variance" of the "occuring energy".

Loudness War is neither dynamic range nor SNR nor variance of occuring engery (although it's related). It's mostly about peak limiting.
For better or for worse, there is a strong desire in the industry to make sure that your song is percieved louder when being played back to back with somone else's song on the radio, a streaming service or a record company executive review meeting.
The song is delivered as a fixed point audio file and this format has a maximum sample amplitude that it can represent. To make it "loud" you need to shove as much energy as possible into the file without any one sample exceeding the the allowable max amplitude. In practice this is done the following way

You record and mix in an unconstrained format (floating point).
The last step is mastering, where the mix session is turned into a fixed point file with whatever format is required (sample rate, bit depth, labelling, streaming constraint, etc.)
During mastering so-called "peak limiting" is applied. You raise the overal level of the song and then shave off the few samples that exceed the maximum allowable amplitude. This generates distortion.
Repeat this until you have achieved the maximum "tolerable" distortion and you can't make it louder without adding too much crap.

In essence you create a trade off between overall loudness and "acceptable" distortion. "Loudness war" basically means that this trade off has skewed hard to towards "loud" and "lots of distortion". A good example to analyze is Tom Petty's album "Hypnotic Eye" the CD version is peak limited, but the high res version is not, so you can actually compare directly what damage the peak limiter has done.
This is quite painful to listen to especially since it's mostly unneccessary. The original version is not as loud, but I can simply adjust the volume to compensate. However, the distortion in the peak limited version is always there.
Back to the original question: you want to look at the "RMS to Peak", since that's exactly what the peak limiters is trying to reduce.  Since these days the "peak" in a file is almost always very close to unity (unless it's "mastered for Itunes"), you can simply take the RMS of the file. Here are some typical examples

0 dB: Square wave
-3 dB: Sine wave
-7.8 dB: "Maggie's Farm" by Rage Against the Machine
-10.8 dB: Averabe of the Billboard 2013 top 100
-15 dB: Average of 10000+ "typical" songs including a wide mix of genre's and release years

Many streaming services don't like this variability since it requires the customer to keep "riding" the volume control and therefore prescribe mastering levels (example -14 dB LUFS which is based on an ITU norm). They may remaster the content or ask the labels to do it which the labels may or may not do. As a result, the same song when played from Spotify may be mastered differently than the same one played from Apple Music.


Microphone Array null steering unexpected results after IFFT I am writing an open source microphone array simulation program:SharpEar.
When I give uniform weight to microphone array elements, the result in the angle is a sinc function as expected. There are 41 elements in array, all equally weighted. Room is  20*20 meters.

Red point in the picture below is a 1000Hz sin wave.  The black areas are filtered by the result of beamforming.

Now I am trying to achive null steering. What I understand from null steering with my computer science brain is: Microphone array will work as a spatial filter.  Fft of the weights of the array will be reflected in angle. That is why result of the uniform weighted microphone array is like the picture above. Now lets say I know there is a noice in -30 degree and the sound that I want to focus is at 0 degree(Angle considered to middle of microphone array red point was at around 0 degree ). And array interests in 90 degrees from -45 to 45. Because I don't want to get anything I am adding a null at -30 degree. And I am having a data like:

0.5, 0.5 , ... , 0 , 0 , 0 , 0.5 ....1, 1, 1, 0.5, 0.5 , 0.5 ...
|                    |                  |
-45              -30 degrees         0 degrees

I believe if I take the inverse fft of the data above and use as weights in my microphone array. I should achive a better focussing in angle domain. But my results are terrible.

I want to ask if my understanding of null steering is correct or am I missing anything?

Ps: Please forgive my lacking descriptions since DSP is not my profession, I am network programmer and trying to run this project in my free time.

I don't understand what you did but I think you didn't understand the idea correctly.
Lets say the weighting vector at a certain frequency is given by $\mathbf{w}^H$ where ${}^H$ is the complex conjugate. That is the elements of $\mathbf{w}^H$ are the weights you apply to each sensor without null steering. And lets say that the array response to the interference is $\mathbf{v}_I$ (which is basically the conjugate of the array weights if you were to point the beamformer to the interference). The simple null steering solves the following optimisation problem:
$$   \underset{\mathbf{w}_{d}}{min}\:\left\Vert \mathbf{w}_{d}-\mathbf{w}\right\Vert ^{2}\\
  subject\; to\\
  \qquad\mathbf{w}_{d}^{H}\mathbf{v}_{I}=0
$$

Using Lagrange multipliers it can be shown that
$$ \mathbf{w}_{d}^{H}=\mathbf{w}^{H}\cdot\left(\mathbf{I}_{N}-\frac{1}{N}\mathbf{v}_{I}\mathbf{v}_{I}^{H}\right)$$
where $\mathbf{I}_{N}$ is the identity matrix and $N$ is the number of sensors.

This solution is not very robust to errors in the interference direction.


Hyperbolic sine "sinh" in fixed-point? As part of an audio peaking filter calculation, I need to implement hyperbolic sine ("sinh") function in a fixed-point arithmetic (fixed-point DSP processor).

What is the proper way to do this?

Be careful with that version of the Cookbook.  I did not actually write it, although I gave permission to Doug to write it.  He has a few typographic errors.  The original has been moved.  I think this is maybe where it lives now but there appears to be a CR/LF problem with the file.

Anyway, to do trancendentals with a fixed-point processor is difficult but possible.  First you need a mathematical description of the transcendental function in terms of arithmetic operations that your fixed-point processor knows how to do.  These are usually an infinite series that is truncated to a finite length.  I have some expressions of these common transcendentals here.

The other thing you have to worry about is scaling of the coefficients to be integer values.  Do you know how to do that in fixed point?  That is fundamental.

Then finally the $\sinh(\cdot)$ is expressed in terms of the exponential function which is described in that Wikipedia article you cited.

Conceptually, doing fixed-point math is not too hard, but it is sorta a bitch.  You must worry about both precision and about overflow.


Filtering an image with two gabor filters I have an image of size 160x160 of pixels. I am trying to convolve this image with two Gabor filters of different scales (7 and 9 respectively) but with the same orientation. I red in an article the following phrase :

The maps (the results of the two convolution operations) have the same dimensionality but they are the outputs of different filters.

After convolving my image (conv2) with these two gabor filters (scales 7 and 9), i got the following results :

The 2 maps (one image after convolving with gabor of scale 7 and the other after convolving with gabor of scale 9) have different dimensionality!

So according to the phrase 1) , why my maps have different sizes ? it is normal ? or false ? i need your opinions.

In other words, i filtered my image with a gabor filter of scale 7, then i filtered my image with another gabor with scale 9 . So i got two maps. Is these maps must have the same size or not ?

If you're using just straight 2D convolution, you will get the full output. What I think you want is conv2(A, B, 'same').

To better understand this, 1D convolution between two sequences with length N & M results in a final sequence of length = N + M - 1. So without any truncation your conv2 result will no longer be of size 160x160, but (160+M-1)x(160+M-1).

If you're using Matlab, see 'help conv2', 'help filter2', 'help imfilter' (if you have the image processing toolbox installed). imfilter defaults to 'same'.

I'll add a little more discussion regarding filter sizes when they contain exponential functions.

Given in reality were sampling a function that extends to infinity in all directions, albeit the values become insignificant at some point, the function does extend that far. The implementation you're referencing fixes the limits to 3 standard deviations, because that's good enough for them. I saw another implementation that used 8. If you're scale, or sigma, is very large, you run the possibility of the filter turning out to be larger than the image you'll be filtering against. What do you do then? Because like I stated earlier, straight convolution produces a sequence that is N + M - 1.

You need to grab only the central part of that result. The central part would also be of the same size as the image your filtering. That's what the 'same' flag does.

One other point I'd like to make is that there's nothing stopping you from just creating a filter that is the same size as your image. Instead of using sigma as the size controller, just specify the dimension of the filter. Given convolution is often times done in the frequency domain, you're going to be zero padding and taking a larger FFT anyway. Look for references on fast convolution using the FFT.

Let's say you're image size is 256x256 and you're filter is 32x32. Here's a simple example without using power of 2 FFT's.

N = 256;
M = 32;
K = N + M - 1; % Length of full convolution

x = randn(N, N);
y = ones(M, M);

% Spatial domain convolution (slow)
r = conv2(x, y, 'same');

% Spatial frequency domain convolution (fast)
z = ifft2(fft2(x, K, K).*fft2(y, K, K));
idx = round((K-N)/2+1):round((K+N)/2); % Indices of the central portion
z = z(idx, idx); % Obtaining the central part of the result

subplot(311), imagesc(r), axis('square'), colorbar
subplot(312), imagesc(z), axis('square'), colorbar
subplot(313), imagesc(abs(r-z)), axis('square'), colorbar


Shifted output of input vs output of shifted input As I understand it, shifted output means first getting the output and then shifting it. We feed in $x[n]$, and get our output as $x[2n]$. Then we shift it by $k$ to get $x[2n-k]$. This is the shifted output.

Now, for output of the shifted input, we shift the input to get $x[n-k]$, then feed it as input to get output as $x[2(n-k)]=x[2n-2k]$

However, the text in the image has these answers reversed. Is there something I'm not getting?

This is all about consequently substituting variables.
Shifting by $t_0$ means replacing the variable $t$ by $t-t_0$. If $y(t)=x(2t)$, then $y(t-t_0)=x(2(t-t_0))=x(2t-2t_0)$.

If you have a shifted input $x_2(t)=x(t-t_0)$, then the output is $y(t)=x_2(2t)$, i.e., the variable $t$ is multiplied by $2$. Replacing $t$ by $2t$ in $x(t-t_0)$ results in $x(2t-t_0)$.

So the slide is correct.

EDIT: Concerning shifting the output, imagine a triangular input $x(t)$ with its apex at $t=0$. The output is a compressed version of the input, i.e., it is also triangular with its apex at $t=0$. If the output is shifted by $t_0$ the apex is shifted to $t=t_0$. The function $x(2t-t_0)$ has its apex at $t$ satisfying $2t-t_0=0$, i.e., at $t=t_0/2$, so the correct formula for the shifted output is $x(2t-2t_0)$.


Simulating a raised cosine filter I am rewriting a simulation written in Matlab to C. I have a problem with a transmit data. The data is filtered by a raised cosine filter. In Matlab it was implemented by using the following function firrcos:
h = firrcos(20,0.25,0.25);

Does anyone know how I can find the documentation of this function ( matlab code realisation)?
Does anyone have a simualtion such filter in Matlab or in C?
EDIT1
Given the order of the filter N=20, carrier frequency Fc=0,25, transition band tbW=0,25 and sampling frequency 8.
I need to initialise such filter in  Simulink.
how to determine a filter span and rolloff factor knowing only the parameters above?
EDIT2
Why the filter impulse response of the filter in Simulink does look different as its implementation in Matlab?
In matlab the upper limit is 0, in Simulink it is 10
EDIT 3
How to initialise h = firrcos(20,0.25,0.25); in simulink?

MATLAB's documentation doesn't recommend the use of firrcos, instead use rcosdesign as recommend. Type doc firrcos  in your command window, here are the first lines:

firrcos Raised Cosine FIR Filter design.
WARNING: firrcos is not recommended. Use RCOSDESIGN instead.

My guess is that firrcos is actually discontinued in new releases. To better help you, try using rcosdesign as recommended. You can actually get the actual MATLAB code for this function by typing the following lines in your MATLAB's command window: edit rcosdesign (you can also check edit firrcos). This opens a new editor window with the actual code for the function. This can help you.
Note that in your EDIT1, Fc is actually the cutoff (not carrier) frequency if you double-check the documentation, and that with three input arguments as you have the default $F_s = 2$. The roll-off factor is a measure of the excess bandwidth expressed as a fraction. The span of the filter in symbols is for you to define.
To add to what has been mentioned by @MBaz, what you need at the transmitter is the root-raised cosine filter. This is so that at the receiver side you also put a  root-raised cosine and the two together result in a raised-cosine filter which is a Nyquist filter with zero ISI while the root-raised cosine doesn't exhibit zeros (in amplitude) at symbol times.
EDIT:
To see the equivalence between the two functions you'll have to give the correct input parameters to both functions. As discussed, first things is having a root-raised cosine. You can do this by setting the 6th and 4th input of firrcos and rcosdesign respectively to sqrt. The equivalents of the remaining parameters are as follows:
beta = R;
sps = Fs/(2*Fc);
span = N / sps;

Then you have the unnormalized two filters as:
h1 = firrcos(N, Fc, R, Fs, 'rolloff', 'sqrt');
h2 = rcosdesign(R, N/(Fs/Fc/2), Fs/Fc/2, 'sqrt');

More on this together with the needed normalization can be found in MATLAB's Communications Toolbox™ Release Notes for R2020b under the chapter on R2013b (I guess it got discontinued after then) titled Functionality Being Changed or Removed.
For your Simulink, click on Simulink Library Browser then on the list on your left choose Communications Toolbox then pick Comm Filters. There you have all your filters, you can set the chosen raised-cosine filter as discussed up.


Effects of averaged periodogram (Welch's Method) on resolution? I'm trying to make an app(running on an android phone) that can detect the frequency of a guitar string being played. I know it'll be hard to detect that but that's a problem I'll deal with later.

Right now I'm trying to reduce the effects of noise on my FFT of an audio signal. I'm sampling at $44.1\textrm{ kHz}$ and have a sample size of $4096$. This gives me a bin size of $10.7\textrm{ Hz}$.

I've read good things about averaged periodograms killing noise. I want to try to use an averaged periodogram, Welch's Method, to reduce the effects of noise.

What I'm wondering is, if I window it into $1024$ sample windows, then will the resolution of my final product become $\approx 40\textrm{ Hz}$?
If I have overlapping windows to keep my resolution high, will I still have effective noise reduction?

Sorry if this is a bad question. This is my first time doing DSP outside MATLAB so I'm struggling.

Yes, if you reduce the FFT window from 4096 to 1024 samples, then you will have 4 times less frequency resolution. Also, if you have significant overlap between the windows, your resolution will not change. It only depends on the length of the FFT window. Furthermore, with Welch's method, a window for each segment is usually used, which again reduces the spectral resolution due to convolution of the signal spectrum with the window spectrum.

So, yes Welch's method reduces noise, but also spectral resolution.

If you want to have both, why dont you increase the recorded signal length to lets say 4*4096? This would amount to 360ms, which, depending on your requirements, might still be fine for guitar pitch detection. THere's always the tradeoff between spectral resolution and required amount of time.


What is the functional difference between DSB and 2SB mixers, and why exactly would 2SB be better for radio astronomy? The Astronomy SE question 2SB mixer in radio astronomy? asked in 2014 in its entirety:

In literature about the instrumentation of radio astronomy, I frequently come across a type of mixer described as "2SB".  I'm familiar with single-sideband (SSB) and double-sideband (DSB) mixers, but I have never studied a 2SB mixer.
What is a 2SB mixer and how does it compare to a DSB mixer?  Ideally, I would also like to know of formal references that directly address the topic for further reading.

is interesting, but its lone 2014 link-only answer in its entirety:

2SB is Dual Sideband, as opposed to DSB - Double Sideband. Here are a couple of papers you might find relevant and interesting:
a 2SB upgrade replacing a DSB1
a paper that mentions the advantages of 2SB over DSB2

leaves me still not really understanding the difference between 2SB and DSB mixers nor exactly when one would choose one over the other.
My impression is that DSB mixers suppress the carrier while allowing both sidebands to pass through a single output, while 2SB mixers further separate each sideband to a separate output.
Question(s):

Is the separation of the two SBs to two outputs the functional difference between 2SB and DSB? If not, what would it be?
The abstract of the second linked paper says that 2SB mixers' "main advantages being to avoid spectral confusion and to diminish effective system temperature by a factor two with respect to double sideband (DSB) receivers." Is it possible explain in fairly simple terms what that is so?

1Hesper et al. A Sideband-Separating Mixer Upgrade for ALMA Band 9 from the 20th International Symposium on Space Terahertz Technology, Charlottesville, 20-22 April 2009
2A calibrated digital sideband separating spectrometer for radio astronomy applications
2Finger et al. A Calibrated Digital Sideband Separating Spectrometer for Radio Astronomy Applications also available at Pub. Ast. Soc. Pac.

Is the separation of the two SBs to two outputs the functional difference between 2SB and DSB? If not, what would it be?

Almost certainly yes, although it seems that the only hits that "2SB mixer" gets is for astronomy sites.  See the figure, from the first paper you cite.  It clearly shows two sidebands being received independently.
Basically, you're getting two SSB channels out of this, which, presumably, can be selected downstream.

The abstract of the second linked paper says that 2SB mixers' "main advantages being to avoid spectral confusion and to diminish effective system temperature by a factor two with respect to double sideband (DSB) receivers." Is it possible explain in fairly simple terms what that is so?

Spectral confusion: DSB mixers take the sideband above the carrier and the sideband below the carrier and mix them equally to baseband.  The output of a DSB mixer is the sum of the downconverted USB and the downconverted LSB.  For uncorrelated signals, the power in the USB and the power in the LSB are added.
Noise temperature: Noise temperature is a function of the ratio of noise power to desired signal power.  Because the DSB mixer output power is the sum of the LSB and USB power, it's output has the noise contribution from both of these sidebands.  Because (based on the context of the article) the desired signal is in only one sideband, you're getting the noise power from two sidebands, and signal from one.
With a SSB mixer (remember you get two SSB channels out of the 2SB mixer) you get the noise power from one sideband, and the signal from one sideband.  So the power/noise ratio is doubled -- this results in the noise temperature being halved.


Acoustic waves propagation distance relation with carrier frequency? In electromagnetic waves, the energy formula is $E=hf$, that means frequency increasing will increase energy and so on the propagation distance.
In underwater acoustic waves, low frequencies can travel larger distances comparing to higher frequencies.

The question is: isn't the relation $E=hf$ is applicable in underwater acoustics?
If it is applicable, does it mean that an increase of carrier frequency will increase the transmission distance of underwater acoustic waves?

For a sine wave of amplitude $A$ and frequency $f=1/T$, its energy from $-T/2$ to $T/2$ is $A^2T/2$. Thus the power depends only on amplitude. The energy over one period does depend on frequency, but in the same observation duration, two sines having same $A$ should have the same energy (assuming observation duration is multiple of their periods).

In communications, we don't transfer energy but signal. Then power, not energy, matters.

Propagation loss is characterized grossly by distance (pathloss model) and absorption that the latter explains the choice of low frequency in underwater acoustics communications. For example, frequencies >= 1 MHz are absorbed very quickly whereas the absorption of low frequencies is much weaker. This is physics and the explaination would be complicated.

FYI, you can use this calculator for sea absorption.

And no, $E=hf$ is for photon, I have not met people who uses quantum mechanics for underwater acoustic waves.


Sketching a tri signal with horizontal shift I'm unsure how to sketch a tri signal with horizontal shift.

Eg , let's say we want to sketch tri((t-4)/4). Can I reduce it to tri( (t/4) - 1 )? This would mean sketching a tri(t/4) signal and shifting it to the right by 1. Or should it be a tri(t/4) signal shifted to the right by 4?

I've found out on wikipedia that the rect signal is defined as rect((X-Y)/T) where you shift the center of the signal by varying Y. But I'm not sure if sketching tri applies the same concept.

The peak of the tri signal is at the point where the argument of tri (that thingy inside the parentheses following tri) equals $0$. The
base of the tri signal extends from the point where the argument equals
$-\frac 12$ to the point where the argument equals $+\frac 12$ (or $-1$ to $+1$ depending on what you are calling tri).  Thus,
tri$\left(\frac{t-t_0}{T}\right)$ is an isoceles triangle whose apex
is at the point $(t_0,1)$ and whose base extends from $t_0-\frac{T}{2}$
to $t_0+\frac{T}{2}$ (or from $t_0-T$ to $t_0+T$ for the other
interpretation of tri).


is it possible to matching b / w images without feature extraction steps? if possible which algorithm may be more suitable? these are example of my image

i have a 245 images similar to it ,i obtained these b/w images after preprocessing stage , i want to use them for recognition persons(not verification) in matlab. can i performe matching algorithm on these image after division into 2 group(one for training and the second for testing) without implementing feature extraction algorithm such as PCA,ICA,.....?

if it is possible , which are the suitable matching algorithms that are preferable for this purpose , especially i need to use matlab ?

any help please?really i need your opinions and consultations.

thanks in advance

Yes, it is possible!

For example, you can take pixel by pixel comparison and apply MSE (mean square error). If the error is close to zero, it is the same picture.

However, you will be able to recognize only those images which are picture perfect matching. For example, you won't be able to account for brightness variations, scale variations, minor transformations, compression loss (e.g. JPEG) etc. So while human being might just see them as SAME image - they are not on every pixel basis. Depending on how much strict or tolerant your application needs to be, do we really need to then look for advanced algorithms that are smart enough to be invariant to these factors.

EDIT:

After your clarification, it does sound like you have more of a special case of matching rather than general purpose recognition problem.

Given the following assumptions:

A gray scale (or color) images are first converted to binary through some thresholding process.
Mostly these images are lines (thick or thin) and the key differentiation is how far these edges are from others.
Edges could be of arbitrary shapes and contours -no geometrical assumptions.
Given the edge detection limitations (in the presence of noise) and thresholding algorithms, edges might be disjoint

Here is what i would suggest to apply:

Edge following algorithms - to make continuous edges where applicable.
Ref #1. http://www.slideshare.net/kiara1011000/edge-following-algorithm-chiara-galdi
Apply Housdourf distances and similar shape/contour matching
Ref #1. Haudorff Distance by Normad
Ref #2. CV Online: The Hausdorff Distance
Ref #3. Hausorff based Matching

Basic idea is that if all edges are more or less at the same locations and of same length, then the images are same.

Given that large portion of region which contains Black and implies no information- MSE would be an extremely BAD idea.

I am only giving you some basic direction toward formulating your problem. As you dig deeper  you can put more specific questions. I would be keen to know what works for you finally.

EDIT 2: Adding one specific solution

One simple possible metric i can think of is:

Say, Ref[i][j] is the reference image and Test[i][j] is the test image.

pixel_similarity = 0;
For each pixel in Ref[i][j]
{
  if ( Ref[i][j] == Edge_pixel )
     (k,l) <- find_nearest_edge_pixel(Ref, Test, i, j, window_size)
     pixel_similarity += weighted_similarity (i, j, k, l)
}

image_similarity_metric = error / (total_edge_pixels_in_ref);

Please Note that this algorithm will be quite limited for heavy scaling and rotation. This is first step. Notify me for any typo or lack of clarity!

Observations:

If you keep two identical images (pixel by pixel)- you will find that (k,l) will be same as (i,j) and weighted error should be = 1.
By dividing the number of edge pixels by total edge population, you will get similarity metric = 1.0 for perfectly matching images.
When edges are almost similar, but slightly misaligned score will decline gradually. The weighted_similarity function should be chosen smartly to deal with what you consider as reasonably aligned edge vs. one that should be discarded.
Other extreme is that if total number of edge pixels in Ref is zero - there is no way to compare it with any other one (since we are considering background pixels having no information related to match).

Try to work this out and see what results you get.

Happy to help you!


MATLAB - How do I get the Filter Function to get the output Energy Spectrum? The image was from a question I posted on reddit, so ignore the title.

I'm trying to write a program that plots these graphs

I know that the first plot is just a plot of the source-signal so the "time" domain of the voice signal. I have that down.
the 2nd plot is the frequency domain signal, which I should be able to get with fft() fcn in matlab? but the plot it gave doesn't look right.

weird fft plot of source signal - me - a female just saying "aaaaahhhh" for about 15 seconds

And after I get the right filter function, I do conv??

Thanks in advance.

What u plotted is the complex fft samples.  In order to see the psd in matlab use plot([-1/2:1/ffteln:1/2-1/fftlen]*fs, 10*log10( abs( fftout)./length(fftout)).^2)).  I am assuming of coarse that your signal is real and therefore plotting from -fs/2:fs/2 is standard practice.  For complex signals usually you want to plot from 0:Fs because nyquist is Fs for complex.


Getting a single feedback peak from a basic 4 stages phaser I try to set a very basic 4 stages phaser (two 2-pole all-pass filters), with a particular feedback path where only one peak is raised at center frequency when the feedback is above 0.

The following amplitude responses describe what I try to get: in the first measure (light green) we get two notches (no feedback). In the second measure (purple) we get two notches and a single peak at center frequency (feedback).

My chain is very simple: the input is sent into two 2-pole all-pass filters having the same center frequency, the output is then mixed 50/50 with the input. I've also added a feedback path, where the output of the last all-pass filter is feed back to the input of the first one (input + state * feedbackAmount).

Instead of getting a single peak I get two (as described in the following picture). I've also tested negative feedback it does not help. Any idea?

Edit1:
My signal flow leading to two peaks when the feedback is close to 1:

Edit2:
My implementation may be wrong indeed, to be sure I've also run a test in Ableton Live 8. Noise straight to the phaser.

First: no feedback:

Second: full feedback, we get two peaks too. (Most phaser I've tested provide the same result):

If I follow the OP’s description properly, the description is what would be a comb filter, given the all-pass filters operate as delay lines (over a certain frequency range of operation). The sum of a signal with a delayed copy of the signal would be a "comb" response, that would cyclically go between constructively summing and cancelling since the phase through a delay line is linearly proportional to frequency (so cylces between 0° and 180° as you sweep frequency). To eliminate the second peak, the bandwidth of the delay line or the feedforward path from the input that adds with the output would need to be limited.

To achieve what is wanted, consider a simple IIR filter with pole placement at the resonant frequency, given by :

$$H(z) = \frac{k}{(z-\alpha e^{j\omega})(z-\alpha e^{-j\omega})}= \frac{k}{z^2-z2\alpha\cos(\omega)+\alpha^2}$$

Where $\omega$ is the normalized radian frequency of the resonance and $\alpha$ sets the bandwidth of the resonance, and $K$ sets the gain. Use $\alpha<1$, with the closer it is  to 1 the tighter the resonance.

For example, with $\alpha = 0.99$, $k=0.1$ and $\omega = 1$ the response is:


Downsampling, Filtering and Averaging I have a signal of bandwidth DC to ,say, 100 Hz. I am interested in the range DC-10 Hz. I oversample by a factor of, say, 30. Thus the sampling rate is 600 samples/sec. Then I bandpass filter and decimate it over a number of stages, say 3 stages (e.g. by a factors of 5, 3 and 2). Thus I get the original signal.

This gives me the advantage of having simpler anti-aliasing filter.

However, there is another advanage of oversampling, which is reducing the noise by taking an average of adjacent (extra) samples. For example, if we oversample by 4, then we take average of every 4 samples.

My questions: Do I need to do filtering before the averaging, as is the case with the downsampling?

EDITED TO CLARIFY.

Averaging is filtering with a rectangular filter kernel, which has a poor frequency response, both in pass band ripple and stop band attenuation, Because of the poor stop band response it won't anti-alias very well.  The non-flat passband will distort the final result even if there is no high-frequency spectral content to alias.


Why is there a dip in this Chebyshev filter response? I want to build and simulate a 6 GHz, 4th order lumped-element Chebyshev bandpass filter with a maximum passband ripple of 0.1 dB and a bandwidth of 1 GHz, and as I understand it, there should be 3 lobes in the S11 response, but instead my circuit has a huge drop at the center frequency. I'm not certain what could be causing this or if this is even common, but I'm assuming it has to do with my implementation. The equations I used to get the component values from the normalized lowpass filter coefficients are as follows and the schematic from MWO is also attached;

The elements in the circuit are (from left to right); [C1 = 3.53078138643953e-06 uF, L1 = 0.199281477338676 nH, L2 = 10.3941804502207 nH, C2 = 6.76935843301265e-08 uF, C3 = 5.63670104271011e-06 uF, L3 = 0.124828215212789 nH, L4 = 6.51082585059285 nH, C4 = 1.08069137002873e-07 uF].If anyone can see where I'm going wrong, and would be kind enough to point it out to me that would be greatly appreciated.

You are confused about which parameters are what.  The $S_{11}$ parameter is the power that's reflected back from the network; the $S_{21}$ parameter is the desirable bit that's passed forward.
Just by conservation of energy, if the filter is lossless then $\left | S_{11} \right | = 1 -  \left | S_{21} \right |$.  So that horrible deep null in $S_{11}$ actually just means that $S_{21}$ is at unity, which is desired.


What is the z-domain transfer function of variable delay? Is it possible to obtain the transfer function of a signal that has a variable delay?

For instance, I have a signal consisting of 10 samples:
0,1,1,1,1,0,0,0,0,1

If I would apply a unit delay to the signal, the transfer function would be z^-1 and the signal is now:

0,0,1,1,1,1,0,0,0,0

What would be the effect e.g. transfer function if I delay the rising edge of the pulse by 1 delay as above, but the falling edge would be delayed by 2 samples? The signal then would be:

0,0,1,1,1,1,1,0,0,0

Could this be derrived as a combination of z^(-1) and z^(-2) or this transfer function is non-linear?

Thank you for your answers.

A (discrete-time) system with a variable delay must be time-varying, and for this reason it cannot be described by a transfer function of the form

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}\tag{1}$$

where $h[n]$ is the impulse response of the system.
Only linear time-invariant (LTI) systems can be described by a transfer function as given in $(1)$.


What is a Phase Delay filter? I'm looking at a job spec, detailed below.   I'm perplexed as to what is being asked for.

Could someone offer a basic understanding of this material?   I can imagine a solution may involve performing an FFT, setting bins outside the desired frequency range to 0, modifying phase of remaining bins, and recombining with iFFT.

but I am left scratching my head as to the purpose of this...

Application is required to:

A) input audio from external microphone jack input,
B) bandpass filter to limit audio to between 50 to 500 Hz
C) add phase delay between 150 degrees to 210 degrees (programmable via settings)
D) output resultant audio stream via speaker phone or audio jack
E) audio output volume to be calibrated to match audio input volume

User interface:
A) two modes- configuration mode, normal mode.
B) in configuration mode allow user to adjust phase delay and other audio    parameters, such as filter characteristics, gain and output volume.
C) in configuration mode user can select audio output via speaker phone or external speakers.
D) ability for user to disable incoming phone calls while app is running.

The thing you're looking for is known as an "all-pass filter", which has a flat amplitude response, but modifies the phase relationships within a signal. This would be a separate stage from the one that limits the bandwidth.

Note that the 50-500Hz bandpass filter is going to modify phase relationships, too, so the all-pass filter will need to be designed so that the overall system response is what the client is looking for.

However, it's possible that the client just wants a time delay, but expressed it in terms of the phase angle of some fundamental frequency. This will need to be clarified.


validation of a pseudo jpeg compression I've done something which is not exactly the jpeg algorithm. What I've done is that given an image $I$ of $N\times N$ size, I have computed the DCT of individual $8\times 8$ blocks.

Once this has been done, I have imposed that the top leftmost corner is left intact but the other elements be equalled to zero.

Then i compute the inverse-DCT on that block and use it to reconstruct my image.

The algorithm seems to work giving the following results

Given that I've ignored lots of steps in the jpeg algorithm (all the things related to huffman encoding, quantization matrix etc) why is this simple algorithm working?

Thanks.

If you compute a $8\times 8$ 2D-DCT, and keep the top left corner only, you are keeping a quantity that is proportional to the average of each  $8\times 8$ block.  This is the DC component, similar to the $0$ frequency in a Fourier transform. This works like a JPEG coding with a quantization table looking like:

$$\displaystyle \left( \begin{array}{cccc}
1 & \infty & \cdots & \infty\\
\infty  & \infty  & \cdots & \infty\\
\vdots  & \vdots &  \ddots & \infty\\
\infty  & \infty  & \cdots & \infty
\end{array}\right)$$

Performing the inverse DCT with this "kind of" quantization spreads the average value over the $8\times 8$ block. What you are doing is equivalent to performing a moving average with an $8\times 8$ uniform window and then downsampling by  $8$ in both the $x$ and $y$ directions. You compress by replacing each $64$-pixel  block by its average, resulting (effectively) in a $64$ times smaller image, and you "fake" an original size image by extrapolating the average.

However, you cannot claim a $1/64$ compression yet, because the coefficients you keep are not quantized yet, and are stored for instance on $32$-bit floating point words. Before you go into more advanced coding, you are about at a  $1/16 = 1/64\times 32/8$ compression rate.


What's the essential bandwidth of the unit step function? The Fourier spectrum is in the Figure, how to find the essential bandwidth?

there can be no essential bandwidth of the unit step function.  defining such would lead you to self-contradiction.

suppose you have a function $x(t)$ which has Fourier Transform

$$ X(f) = \int\limits_{-\infty}^{+\infty} x(t) e^{-j 2 \pi f t} \ dt $$

and from $X(f)$, you defined some consistent measure of bandwidth $B_x$ in such a way that is independent of the amplitude of $X$ and is dimensionally consistent with $f$.

now suppose you scale the time argument of $x(t)$ to speed it up or slow it down:

$$ y(t) = x(a t) \quad \text{for } a > 0 $$

we know that the the Fourier Transform of $y(t)$ would be

$$ Y(f) = \frac{1}{a} X\left(\frac{f}{a} \right) $$

and the bandwidth of $Y(f)$ defined in the same manner as $B_x$ would be

$$ B_y = a B_x $$

so, suppose $a>1$, then you speed up $x(t)$ by the factor $a$ and the bandwidth is increased by the same factor.

but if you "speed up" (or slow down) the unit step function, you still have the very same unit step function:

$$ u(t) = u(a t) \quad \text{for } a > 0 $$

where

$$ u(t) \triangleq \begin{cases}
1, & \text{if } t > 0 \\
0, & \text{if } t < 0 \\
\end{cases} $$

so scaling time for the unit step does not change it, yet if it had a bandwidth, scaling the time must change that bandwidth.  a contradiction results.


DFT and then IDFT does not provide the same signal I'm new at Fourier Transforms. I'm using AForge.Math FFT/DFT algorithm (C#). Said that I have this situation:

My signal:

$$
1.4 + sin(2\pi t) + 0.4sin(5·2\pi t)
$$

Sampling: 100 samples each second.

Here is the spectrum. The DFT finds the frequencies of my signal: 1Hz and 5Hz (and the constant component: 0Hz).

Here is the inverse transform. As you can see the result is similar in shape to the original but obviously something is happening there.

My questions:

What parameters can I change to get better results.
Should I expect to completely be able to regenerate the signal in time domain from the spectrum?
Can the IDFT be used as a signal generator? It's very easy to create the peaks at the desired frequencies.
Maybe I need to go to the books, but, why there is always a mirrored frequencies in the DFT and how should I interpret them? I don't have any frequencies at 95Hz nor 99Hz. Is it a flaw of the DFT?
Why the initial point is not mirrored?
The amplitudes of the peaks in the spectrum are not correct: the first one should be 1 and the second 0.4. I've heard that the amplitudes in a spectrum are not important, is that truth, if so why? Anyway, I think it's very useful to know what is the amplitude of the sinusoids, is there a way to do that?
Can I throw away the symmetric part?
If so, can I throw away the almost null values till the first frequency from the middle to start? (in my case from 50Hz to almost 5 Hz).

EDIT:

As pointed here I needed to, first keep the real/imaginary data from the direct transform and, second, adapt that data because I'm using an complex inverse DFT.

Here are the results:

In cyan/purple : Re/Im parts

In red: final result.

The X spacing is wrong but that's a story for another day.

Thanks everyone.

If you want your IFFT result to match the input to the FFT, don't throw away either the "imaginary" (or sine correlation) components, or the "upper half" of the spectrum, before the IFFT.  Also make sure the input to the IFFT is exactly conjugate symmetric, so that all the imaginary components in the IFFT results get cancelled out (except for rounding errors or quantization noise).  The length of the IFFT has to be the same as that of the FFT.

Your 99 and 95 Hz components are also the -1 and -5 Hz components needed (in complex conjugate symmetry) to cancel out any rotation of the time domain data into the imaginary plane.  Since the FFT is a complex to complex transform, you can't just ignore these components even when dealing with strictly real data.

The 0 Hz or DC component of an FFT result appears double height because it includes it's own mirrored value.

Unless your sinewave's are exactly periodic within the FFT length, you will need to interpolate the true amplitudes, since their peaks will be between the FFT result bin centers.

The IDFT and IFFT can be used as signal generators, although great care must be taken with phase for any frequencies that are not exactly periodic in the FFTs length.  See Phase Vocoder analysis/synthesis for details.


Double Threshold Mechanism I was studying an article about optical music recognition. Somewhere in segmentation section is mentioned that

  The final step consists of isolating the note heads. This is performed by using a double threshold mechanism to obtain the results shown in Figure.

I couldn't find anything on the internet that explain what does "double threshold mechanism" means. I guess it is a very basic question. I hope someone guides me.

without knowing the context, it could be a number of things but one does not threshold unless you are making a decision (like a hypothesis test).

One example of a double threshold is getting through security at an airport. You go through the metal detector. The machine alerts or doesn’t . That is the first decision. It has a high false alarm rate because the cost of a missed detection is high. It alerts and you empty your pockets (hoping no one is stealing your laptop) and the machine still alerts.  The pat down is the second decision.  The pat down takes longer and has a higher cost but is usually more accurate, unless the gate is run by perverts, but that is something else.

There are usually more thresholds that you need to pass at the airport, but the general idea of using multiple “independent “ tests to achieve better classification performance holds.


frequency translation or pitch shifting through linear interpolation in frequency domain I want to perform pitch shifting with the method described in NEW PHASE-VOCODER TECHNIQUES FOR PITCH-SHIFTING, HARMONIZING AND OTHER EXOTIC EFFECTS. Basically, it involves:

Locate the spectral peaks and their indices in spectrums.
Determine the region of influence (ROI) for each detected peak.
Perform peak shifting according to the amount of frequency shift. Specifically, peaks and their ROI are shifted by an amount of frequency.
Update the phase to maintain the phase coherence, as in phase vocoder for time-scaling.

I am stuck in step 3. The paper describes that with a fractional number of frequency bin shifting (rather than integer number of frequency bin), one should use interpolation since the frequency bins are discrete:

  If the amount of frequency shift    is a fractional number of frequency bins, then frequency-domain interpolation is required, since the sinusoidal peak is only known at discrete frequencies.

My problem is that how do I perform peak shifting through, say, linear interpolation? And, assume a shifted spectrum is obtained, as interpolation insert values between the original discrete bins, how is it finally istf to the pitch-shifted output? I appreciate for any help.

EDIT:

I've realized my question can be simply be phrased as "fractional shifting of spectrum components". The following code is an example to achieve fractional shifting with a random signal using linear interpolation:

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(5)
y = np.random.random(20)
alpha = 3.2    # shifting factor

# Pad the signal with zero
# to avoid spill-out caused by
# large shifting factor
nPad = np.ceil(alpha*len(y))
yPad = np.concatenate(
           (y, np.zeros(int(nPad))))

x = np.arange(len(yPad))
xInt = x[0:len(y)]*alpha    # the 'xp' in function interp
yInt = y                    # the 'fp' in function interp

# Now define 'x' in function interp
start = np.ceil(xInt[0])
end = np.ceil(xInt[-1])
if float(xInt[-1]).is_integer():
    xQue = np.arange(start, end+1, 1).astype(int)
else:
    xQue = np.arange(start, end, 1).astype(int)

yNew = np.interp(xQue, xInt, yInt)    # linear interpolation
shifted = np.zeros(len(yPad))
shifted[xQue] = yNew

plt.plot(yPad, 'b')
plt.plot(shifted, 'r')
plt.show()

The above example uses a randomly generated signal to illustrate fractional shifting of a spectrum.

However, in the real case, a spectrum shouldn't be padded with zero in the first case as I do in the code. So, how to deal with the "spill" with spectrum component if the shifting factor is too large?

You can use a Sinc (or windowed Sinc) kernel for interpolation between FFT/DFT result bins.  Note that near DC and Fs/2, the interpolation needs to be circular.

After shifting a peak up and determining its phase, you can again use Sinc interpolation to figure out how it should be represented in some number of nearby (non-fractional) bins.

Added:

Linear interpolation can lead to the distortion of harmonic ratios.

To avoid aliasing, the result has to be low-pass filtered or bandlimited to below Fs/2.  So any spectrum shifted above that needs to be filtered out and thrown away.


Relationship between the IDFT of a sampled DTFT and its discrete-time domain signal Suppose we are given an input signal s[m,n] with DTFT $S(\omega_1, \omega_2)$.

We sample it at $\omega_1 = \frac{2 \pi k}{256}$ and $\omega_2 = \frac{2 \pi l}{256}$ to get a 256 point DFT S[k,l]. Now suppose we take the IDFT of S[k,l] to get $\tilde{s}[m,n]$

I am trying to understand the relationship between $\tilde{s}[m,n]$ and s[m,n]. My understanding is that I should be able to express $\tilde{s}[m,n]$ as a summation of every 256th sample of the original signal.

So far I have been advised to try expressing the Sampled DTFT as a dirac comb, i.e.:

$S[k,l] = \sum_{k = -\infty}^{\infty} \sum_{l = -\infty}^{\infty} S( \frac{2 \pi k}{256}, \frac{2 \pi l}{256}) \delta(\omega_1 - \frac{2 \pi k}{256}, \omega_2-\frac{2 \pi l}{256}) $

and use the definition of the IDFT:

$\tilde{s}[m,n] = \frac{1}{256^2}\sum_{k=0}^{255} \sum_{l=0}^{255} S[k,l] e^{j(\frac{2 \pi km}{256}+\frac{2 \pi ln}{256})}$

Using these two pieces of information I can get the expression:

$\tilde{s}[m,n] = \frac{1}{256^2}\sum_{k=0}^{255} \sum_{l=0}^{255} [\sum_{k = -\infty}^{\infty} \sum_{l = -\infty}^{\infty} S( \frac{2 \pi k}{256}, \frac{2 \pi l}{256}) \delta(\omega_1 - \frac{2 \pi k}{256}, \omega_2-\frac{2 \pi l}{256})] e^{j(\frac{2 \pi km}{256}+\frac{2 \pi ln}{256})}$

but I am not sure how to proceed from here. Any advice would be helpful I have been racking my brain over this for a couple hours now.

For simplicity, I prefer a 1D notation without losing generality.

The IDFT associated with the properly taken uniform (in frequency) samples of a given valid DTFT can be shown to be the following.

Case-1: Assume a finite length signal $x[n]$ of length $N$, and consider its DTFT $X(\omega)$. Now if you take $M \geq N$ uniform (frequency) samples of $X(\omega)$ according to $$ X[k] = X( \frac{2\pi}{M}k) ~~~, k=0,1,...,M-1 $$ then the M-point inverse DFT of $X[k]$ will exactly be $x[n]$ padded with $M-N$ zeros; i.e.,
$$ y[n] = IDFT_M\{ DFT_M\{ x[n] \} \} = x_M[n] $$ So the sequence $x[n]$ and $y[n]$ have the same values for the first $N$ samples (i.e., all samples of $x[n]$ are retained in $y[n]$), but $y[n]$ have further (effectively padded) zero samples.

If $ M < N$ , then time-aliasing in $x[n]$ will happen and the first $N-M$ samples of $y[n]$ will be in error. But the remaning $2M-N$ samples of $y[n]$ will be identical to the original $x[n]$ at the corresponding index positions.

Case-2: $x[n]$ is infinite length (and stable) sequence and $X(\omega)$ is its DTFT. Then according to $M < N$ subcondition of case-1, all samples of $y[n]$ will be aliased. No samples of $y[n]$ will be equal to $x[n]$, except other than (pure) coincidence.


Difference between frequency sampling and windowing method These two filter design methods seem fairly similar. Both ask to determine a frequency response and then take a inverse DFT to get the filter coefficients.

So what's the difference?

The Windowing method for filter design uses samples of the impulse response for the filter coefficients, truncated to the filter length (and then possibly, but not required, tapered by a higher performing window).

The Frequency Sampling method for filter design uses the IDFT (Inverse Discrete Fourier Transform) of the desired frequency response for the filter coefficients, using the number of samples across the frequency response to be equal to the number of coefficients in the filter.

The point of confusion is that in practice the impulse response is often approximated by an IDFT of the target frequency response, while an IDTFT (Inverse Discrete Time Fourier Transform) would actually be required. The IDTFT can be approximated by interpolating the frequency response first prior to taking the IDFT such that the time span of the impulse response increases. (Allowing the impulse response to sufficiently decay).

Therefore, if you take the IDFT of the sampled frequency response with the number of samples equal to the number of filter coefficients, you are using the Frequency Sampling method. If you take the IDFT of a higher number of samples in the frequency response, and then truncate to your filter length, you are approximating the Windowing method to the extent that the impulse response matches the actual impulse response for your target filter.  (For rectangular response filters there is no need to do an IDFT for the Window method; just sample a sinc function since that it easy to calculate and represents the actual impulse response).

(For further description of how the DFT and the DTFT are related see my response to:
Sine Wave Aliasing during IFFT)

This is a useful reference as well:

http://blogs.mathworks.com/steve/2010/03/15/the-dft-and-the-dtft/

Frequency Sampling

Very simple. Not very effective (except where selective precision is needed).

Coefficients of the FIR filter are the inverse DFT (discrete Fourier Transform) of the sampled frequency response desired.  This approach is very simple, but is not very effective for typical filter designs; much longer filter lengths are required over optimized approaches such as least squared and equi-ripple (Remez-exchange) filter design approaches. To be noted, the solution will be exact at the frequency sample locations, with ripple in between. When an application requires an exact solution at discrete frequency locations (and you do not care what happens in between) this can be a practical approach.

Windowing

Simple. Non-optimal.

Coefficients of the FIR filter are determined by an approach that approximates the inverse DTFT (Discrete-Time Fourier Transform); specifically the frequency response is interpolated (as a true IDTFT would be of a continuous frequency function), an IDFT is actually still calculated and the result of the IDFT (samples of the impulse response) is windowed by multiplying by a Window function of choice (I always liked the Kaiser window with its two "knobs" N and Beta). Multiplying the impulse response by a rectangular window is the estimation of the IDTFT process; as zero padding a time domain function will result in additional samples of the DTFT in the frequency domain.  The frequency domain is interpolated such that the impulse response determined from the DFT decays to insignificant values. Note that this results in the time domain samples that are samples of the actual ideal impulse response, simply truncated in time (this highlights the second view of what the difference is between the two approaches; here we see that prior to truncation and windowing, the impulse response used is samples of the continuous time impulse response since we over-sampled in the frequency domain to eliminate the effects of time domain aliasing.

Demonstration Example

Here is a specific example demonstrating the difference between Frequency Sampling and Windowing approach. Consider a target filter design with 99 taps and the following desired response for the normalized frequency range of $[0: N-1] * 2\pi/N$:

desired_resp= [ones(1,10) zeros(1,80) ones(1,9)];

For the frequency sampling method, the coefficients of the filter are simply the IFFT of the desired resp (rotated such that the dominant tap is centered in the filter, which would have naturally occurred in the IFFT if we had included the expected filter delay).

As shown in the frequency response for these coefficients, this will result in an exact solution at the frequency sampling locations (which can be useful in certain applications), but due to the fact that the impulse response has not fully decayed at the edges of the filter, there will be significant ripple in between. This also means that the impulse response shown is not truly samples of the impulse response for the target frequency response, as we are seeing the effects of folding in the time domain (equivalently to the effects of aliasing in the frequency domain when sampling in the time domain- it is the same process). This is specifically what differentiates the Windowing method of filter design!

In contrast, the Windowing method uses 99 samples of the actual impulse response. We can approximate this by interpolating our desired frequency response (at a higher sampling rate in frequency, which corresponds to a longer time span in time), or certainly in an exact approach, we can mathematically calculate the impulse response and sample that directly, skipping the IDFT operation all together [That in fact would be the exact definition of the Window approach, but in practice it can be sufficiently approximated by a IDFT if the target response is sufficiently interpolated. To demonstrate the difference, I have interpolated our target response by a factor of 10, resulting in the following estimate of the impulse response (989 samples) by taking the IDFT of the desired frequency response after it has been interpolated (note the result still suffers from time domain aliasing compared to using the formula for a sinc function as the impulse response which would be precise, but this is to show that even though both solutions use the IDFT, the result of taking an IDFT of 99 samples directly versus taking 99 samples of a 989 sample IDFT are very different):

Below is the comparison of the coefficients in both approaches showing how the difference in values occur at the tails; the Window method better approximates the actual impulse response, which is still decaying past the length of our filter. The Frequency Sampling method causes the values beyond the length of the filter to fold over, increasing the values of the coefficients, leading to increased ripple in the frequency domain.  Thus truncating the result (rectangular Window) using the Window method results in a higher performance filter, and further reducing the tails with a more advanced window function can lead to further performance at the expense of increasing the transition band.

The center 99 samples of this impulse response are chosen as the filter coefficients (this is the rectangular window). Below shows the result for both filter designs. Further, the performance of the windowing approach can be further improved by selecting the 99 taps used with a higher performance window such as Hann, Kaiser, Blackman-harris, etc.

Zooming in on the passband further shows how the frequency sampling method results in an exact answer at the frequency sampling locations but has overall more error than the window approach with the same number of taps.

Note on FIR Filter Design related to the Windowing and Frequency Sampling Approach

I believe the above two approaches are insightful, but I would rarely use either for an actual filter design. It would add interest if readers commented on applications where either of these two approaches are actually the "best choice". My design approaches of choice are the Least Squared and the Remez-Exchange (equi-ripple) algorithms. Most of the time I go with Least Squared as that results in the best statistical performance (best SNR solution), while the Remez will meet a target specification consistently throughout the passband and stopband.

See below plots of comparative results for a 31 tap FIR filter solution with a 9.6KHz passband and 12KHz stopband showing how a Window design (using the Hamming Window) compares to Least Squared, and how Equi-ripple (Remez) compares to Least Squared; giving insight into why I typically go with Least Squared.

Also see Matt L's response to Difference between Equiripple & Least Squares design for FIR digital filters.


Reduce the Number of Intensity Levels of a Grayscale Image in MATLAB I have written a Matlab script to reduce the number of intensity levels of each pixel of a grayscale image from 256 to some power of 2.

img_color = imread('photo.jpg');
img_gray = rgb2gray(img_color);
imshow(img_gray);
[rows, cols] = size(img_gray);
noOfDesiredIntensityLevels = 2; // test data. will check for 4,8,16,32,etc.
bitsNeededToRepresentIntensityLevels = log2(noOfDesiredIntensityLevels);
new_img = img_gray;
for i = 1 : rows
    for j = 1 : cols
        new_img(i,j) = floor(img_gray(i,j)/(2^(8-bitsNeededToRepresentIntensityLevels)));
    end
end
figure
imshow(new_img);

On execution, the script returns a black image. My expectation was that the image will be turned into black-and-white (intensity value for each pixel will be either a 0 or a 1).

What am I missing here?

P.S: I am a novice in Matlab and Image Processing. So, please ignore any mistakes in my understanding.

I think by number of levels you want the image full scale grey-scale to be divided piece-wise into given number of levels.
For example: -

If number of levels = 2, then you want only two grey-scales in your image i.e. (0 and 128) or (128 and 255) depending upon if you are using floor or ceil within the range
If number of levels = 4, then you want 4 different levels in your image i.e. (0, 64, 124, 192)

Solution -
This operation can be done in a single line

new_img = ceil(img_gray./step)*step;

Here I have divide the entire range into desired number of parts (levels) with size as step and used ceil function to restrict the result to the upper bounds (you can use floor for lower bounds)

Full Code -

img_color = imread('peppers.png');
img_gray = rgb2gray(img_color);
imshow(img_gray);

[rows, cols] = size(img_gray);
noOfDesiredIntensityLevels = 2;
step = ceil(255/(noOfDesiredIntensityLevels - 1));

new_img = ceil(img_gray./step)*step;

figure
imshow(new_img);

% optional code to show the levels
allSteps = 0;
currStep = 0;
while(currStep < 255)
    currStep = currStep + step;
    allSteps = [allSteps currStep];
end

allSteps

Edit 1 - I have included few extra lines which shows you the different levels (steps) in an array. Its just optional but nice to know the exact values of the levels.