Prerequisites for audio signal processing I started a DSP course online and quickly discovered that as much as some people say that it is possible to avoid the more advanced mathematics of signal processing, this didn't seem the case.

It seems like a solid understanding of Calculus and Probability Theory is needed to properly grasp the subject and to be able to be creative with it in the audio field.

Which parts of Calculus are mainly used for dsp? Are there any parts than can be excluded?

DSP is indeed a very mathematical discipline, but not completely. The amount of mathematical knowledge you need to operate DSP concepts comfortably is luckily limited to a particular fairly small subset. I would say that you need the following:

Complex variables. DSP deals with oscillating signals and systems, and those are very conveniently represented using Euler's complex exponential formulas. This includes knowledge of polynomials over complex variables, and so on.
Simple differential calculus. I would say that you need to know how to take derivatives of polynomials and other common functions such as exponentials and logarithms. Knowing l'Hôpital's rule also occasionally comes in handy. Taking derivatives is primarily required for finding minima and maxima of functions. Note that you need to be able to do this for complex variables as well, but at this stage there is practically no difference in concepts between real and complex differentiation.
Linear algebra. I put linear algebra before integration because it is much more useful, and since integration is a linear operation, a lot of it can be replaced by linear algebra concepts. I can't exaggerate how important linear algebra is, so go and learn it. I highly recommend Gilbert's MIT Open Courseware linear algebra class. No one teaches it better than this guy. Just watch the lectures.
Integral calculus. I'm not quite sure how to approach describing the importance of integration in DSP. If linear algebra was a moderately easy, fun and useful topic to learn, integration, especially complex integration, takes a lot of commitment and intuition, both of which require lots of time. Even though I usually advocate for in-depth mathematical knowledge, here's what you can very easily get away with. Learn about integration concepts, what it means and how integration is related to derivatives. Since most integrals are difficult to work out by hand, people use integration tables to look up answers instead. This will get you through most integration you need to know for DSP.
Infinite series. This is the one part of calculus that everyone forgets the minute they take that final exam, but it actually comes in handy in DSP. Specifically, when analyzing recursive systems, i.e. those that take their output and reuse it as input creating a feedback loop, you need to be able to spot several types of infinite number patterns that appear and point to a single formula that produces every single term in that infinite string of numbers or polynomial coefficients.
Probability. I hesitate a little to write about this one, because not all DSP courses strictly require probability, but if the one you're looking at does, it wouldn't hurt to know it. Not be be lost, you need to know expressions and concepts of means and variances, as well as Bayes' rule.

Obviously, what I've given above are the minimal requirements, but you may need to learn more as you go.


chirp spread spectrum with 2 users I am a beginner in digatal communications. I have 2 users and I want use a chirp spread spectrum.

Can we take an up-chirp and a down-chirp as spreading signals for the 2 users? and if yes, why? The down-chirp is $x_{down}(t)=x_{up}(-t)$ with $x_{up}(t)$ the up-chirp.

Yes, you could. When you're using spread spectrum for multiple access, what's important is the cross-correlation between any two users' spreading signals. Any cross-correlation between two users' spreading sequences will cause "leakage" from one user's channel into another's; this causes an artficial noise floor that will limit your system's performance.

You'll find that an up-chirp and down-chirp have very little correlation with one another, so they would be usable for this purpose.


Why cosine is not an eigen signal? According to this website:

  If the output of a system has the same type as its input signal, then the input signal is referred to as the eigen function of the system.

but in this question it is saying only option $(c)$ is eigen function.

My question is why can't it  be option $(b)$ and $(d)$ as it will also obey the defintion of eigen function as given in the website?

PS: Please let me if know my question is not clear i will edit it but please don't downvote.

That definition is rubbish.

An eigenfunction of a system is a function that passes that system, only scaled by a scalar value.

Now, suppose your system is a simple complex band pass filter, which has a pass band at $f_0$, but stop bands everywhere else, especially at $-f_0$.

Now your cosine $\cos(2\pi f_0 t)$ has spectral components at both $f_0$ and $-f_0$. The output of that filter contains only a multiple if the positive part. Hence, the output is not a scalar multiple of the input cosine. Done.


Why do we need to minimize the TX power to require for meeting QoS? Here is the part of the paper,

  Based on these architectures, a lot of recent developments have lately appeared intending at enhancing the rate-energy performance of multiuser Multiple
  Input Single Output (MISO) SWIPT systems.These works mainly target at the optimization of TX precoding and PS operation, and can be classiﬁed into the following two categories. The ﬁrst category is based on whether RXs are required to perform both Information Decoding (ID) and Energy Harvesting(EH) (co-located ID and EH) or just act as ID or EH RXs (separated ID and EH). The second category includes performance objectives like the minimization of TX power required for meeting Quality of Service (QoS) and EH constraints, and throughput or EH maximization for a given TX power budget and QoS constraints

Why do we need to minimize the TX power to require for meeting QoS? Does't the more powerful signal means it can be broadcasted to the farer place?

paper link : https://www.groundai.com/project/energy-sustainable-iot-with-individual-qos-constraints-through-miso-swipt-multicasting/

I interpret that sentence the following:

  Under energy harvesting contsraints and Quality of Service constraints, the second category of works strive to minimize TX power.

So, TX power isn't minimized in order to meet QoS constraints, but as far as possible while still meeting QoS constraints.


What is the peak resonance of convolving with a sine FIR filter? I'm trying to improve my understanding of FIR filters. As an experiment, I've manually created an FIR filter, whose coefficients follow exactly one period of a sine wave. I'm wondering what is the peak gain frequency of such a filter.

For example, I'm filtering a 44.1 kHz audio signal with filter coefficients corresponding to a 440 Hz sine wave1. Intuitively I would expect that the convolution reinforces this frequency. However, I noticed that the frequencies of slightly lower frequencies (roughly a few semitones, ~365 Hz) actually have a stronger gain. This is also confirmed by running the filter coefficients through scipy.signal.freqz (the red line corresponds to 440 Hz, the gray lines indicate frequencies of neighboring semitones):

Most likely I've encountered a basic rule of digital filter design that I'm not aware of. I'm wondering:

Why is the peak shifted w.r.t. the frequency that is used in the convolution?
Is there a formula that describes the relationship between the two, i.e., is it possible to compute the frequency that the sine convolution should have so that the frequency response peaks at 440 Hz?

1 With a sampling frequency of 44.1 kHz, the period of a 440 Hz sine has a fractional length 100.23 samples, which I'm rounding down to 100 filter coefficients. Could this imperfection play a role in the effect I'm seeing?

If your signal was a windowed complex exponential, then the maximum of the magnitude of its discrete-time Fourier transform (DTFT) would equal the frequency of the complex exponential:

$$x[n]=\begin{cases}e^{jn\omega_0},\quad&n\in[0,N-1]\\0,\quad& \textrm{otherwise}\end{cases}\tag{1}$$

The DTFT of $(1)$ is easily computed from the convolution of the DTFT of a rectangular window $w[n]=1$ for $n\in[0,N-1]$ and zero otherwise, and the DTFT of $e^{jn\omega_0}$:

$$\textrm{DTFT}\{w[n]\}=W(e^{j\omega})=e^{-j(N-1)\omega/2}\frac{\sin\left(N\omega/2\right)}{\sin(\omega/2)}\tag{2}$$

$$\textrm{DTFT}\{e^{jn\omega_0}\}=2\pi\delta(\omega-\omega_0),\qquad |\omega|<\pi,\quad |\omega_0|<\pi\tag{3}$$

Note that the DTFT in $(3)$ is of course $2\pi$-periodic.

With $(2)$ and $(3)$ we can write the DTFT of $(1)$ as follows:

$$\begin{align}X(e^{j\omega})&=W(e^{j\omega})\star\delta(\omega-\omega_0)\\&=W\left(e^{j(\omega-\omega_0)}\right)\\&=e^{-j(N-1)(\omega-\omega_0)/2}\frac{\sin\left(N(\omega-\omega_0)/2\right)}{\sin((\omega-\omega_0)/2)}\tag{4}\end{align}$$

The magnitude of $(4)$ clearly has a maximum at $\omega=\omega_0$ $(+2\pi k)$

However, for a windowed sinusoidal signal we get two contributions of the form $(2)$ centered at $\pm\omega_0$ because $\sin(n\omega_0)=(e^{jn\omega_0}-e^{-jn\omega_0})/2j$:

$$y[n]=\begin{cases}\sin(n\omega_0),\quad&n\in[0,N-1]\\0,\quad& \textrm{otherwise}\end{cases}\tag{5}$$

$$Y(e^{j\omega})=\frac{1}{2j}\left[e^{-j(N-1)(\omega-\omega_0)/2}\frac{\sin\left(N(\omega-\omega_0)/2\right)}{\sin((\omega-\omega_0)/2)}-\\\quad e^{-j(N-1)(\omega+\omega_0)/2}\frac{\sin\left(N(\omega+\omega_0)/2\right)}{\sin((\omega+\omega_0)/2)}\right]\tag{6}$$

I don't think that there is an analytical formula for the location of the maximum of the magnitude of $(6)$, but it is clear that due to the combination of the two terms in $(6)$, the maximum shifts away from $\omega_0$. That shift can become especially large for frequencies either close to DC or close to Nyquist. It is smallest for frequencies around half the Nyquist frequency (i.e., $f_s/4$). For small frequencies $\omega_0$, the actual frequency where the maximum occurs is smaller than $\omega_0$, whereas for large frequencies it is greater than $\omega_0$.


Affine Transformation for matching 2D point sets I have two sets of 18 corresponding 2D points. These coordinates represent the main joint points of a human body, describing the "human pose".

I want to decide if these two poses are similar/have the same shape. In order to do so, I find the affine transformation (matrix) that transforms the second pose to the first pose (using least-squares). This gives me the image of the second pose pictured on the first one.

Now, my question is, what is the best way to conclude if the image of the second pose is similar to the first one? I'm thinking about using a distance metric like mean square or SAD, but I'm not sure. + What about normalizing/standardization?

There are two ways by which you can approach this problem:

Match the point clouds directly with RANSAC.
Consider matching the resulting graphs (closely related to 1, still)

What you have done already is closest to #1. To complete it, you need to repeatedly obtain subsamples, find the transformation and then check its "error". The "error" is how many projected points are within a distance threshold (which you are going to have to estimate). For more information, please see this thesis. It does a great job at defining the problem in the beginning (which is very important) and then talking about the approaches.

But, since you are already dealing with pose estimation, you can incorporate its constraints into the model and go for estimation via graph matching. In this case, you would be looking at extracting a graph from your images and then compare it with the pose graph, OR, compare the extent to which two graphs coincide (the graph obtained from the left image to the graph obtained from the right image). A key part of this is the derivation of image skeletons. Gonzalez' book has a very good section on skeletons (chapter 11) and their applications.

The existence of a graph imposes certain constraints on the relationships between the points. For example, the hand "leafs" are expected to be no more than 5-7 "hops" from the root node (e.g. the head node) and this can improve the detection rate. For more information please see this paper. Or more generally look for "graph based pose estimation" and "graph recognition". Be aware that "graph matching" as a term is associated with something entirely different in mathematics.

Hope this helps.


Beginner Educational Resources New member here.  Recently got hired into a mid level EE position.  Have a BSEE but performed more radar performance analysis at a higher system level previously.  In my new group we're doing a lot of digital signal processing design, validation and verification.  I have a strong background with Matlab and am looking for educational resources on pulse doppler, pulse compression, a/d converters, FPGA design using verilog.  Appreciate any input that anyone can provide.  Thanks guys and process on.

Any of the Stimson airborne radar books are good modern resources for general radar signal processing.
For more specific applications (basic STAP, EW, mismatched filters etc), you can look at Mark Richard’s book(s).
If your company has a library, check to see what radar books they have on hand; you may be able to get these for free.
Last, MATLAB has entire packages based around radar systems and plenty of examples that use them. If you have access to a MATLAB license, those might be good to check out as well.


Convolving complex exponential with box function (discrete) If I have $$x[n] = u[n] - u[n-4]$$ where $u[n]$ is the unit step function, and $$h[n] = e^{-i\frac{\pi}{2}n}$$

does $x[n] * h[n] = 0$?

I tried doing the convolution sum and I got: $1 - i - 1 + i = 0$ and I also just tried plotting points and think I got 0 also, but I'm not sure if I did either of these correctly. Is there any easier or intuitive way to think about convolving a discrete complex exponential with a function, or (if it's easier the other way around), convolving a discrete box with another function?

You can use the convolution theorem to convert the convolution to the product of Fourier transforms of the given signals.
$$h[n]*x[n]=\mathcal{F}^{-1}\{X(\omega)H(\omega)\}$$

$$H(\omega)=\mathcal{F}\{h[n]\}=\delta(\omega-\frac{\pi}{2})$$
$$X(\omega)=\mathcal{F}\{x[n]\}=e^{-j\frac{3 \omega}{2}}X_r(\omega)$$ where $$X_r(\omega)=\left\{\begin{matrix}
\frac{\sin(2\omega)}{\sin(\omega/2)} &\omega\neq 0 \\
4 & \omega=0
\end{matrix}\right.$$

However since $X(\omega)\delta(\omega-\frac{\pi}{2})=X(\frac{\pi}{2})\delta(\omega-\frac{\pi}{2})$ and $X(\frac{\pi}{2})=e^{-j\frac{3\pi}{4}}\frac{\sin(2\pi/2)}{\sin(\pi/4)}=0$, the inverse Fourier transform and therefore the covolution is zero.

You can also assume $x[n]$ is the impulse response of an LTI system and $h[n]$ is an input signal. Since the complex exponential is eigenfunction of the system, the output whould be an scaled version of the input i.e. $Ah[n]$ where $A=X(\frac{\pi}{2})$ and $X$ is the FT of $x[n]$.


How to window a time domain signal for 6 symbols with cyclic prefix length of 512 samples? Based on LTE stds. I am running a simulation for 6 symbols (each having 1200 central subcarriers) each in a 2048 FFT/IFFT bin. To these symbols are appended 512 samples as cyclic prefixes. I have to window each symbol to see the impact on OOB. How do i window them ?

Assuming that it would suffice to have an average result of all six symbols, I suggest that you simply concatenate all symbols with cyclic prefix in time and then window and FFT the whole dataset. This would give you the same result as an average of the result you would get for each symbol + prefix independently. Therefore the window length would be (2048+512)*6 = 15,360.

You would only consider overlap techniques when you need to continuously stream an fft process, or if you have a dataset that is not feasible to process with one fft. 15,360 is quite reasonable to do as one fft process.

Windowing involves multiplying the time domain data by the window weight sample for sample. After windowing the resulting windowed data is FFT'd. Purpose is to reduce spectral leakage from strong tones, thereby increasing the dynamic range and ability to see weaker signals together with stronger signals. This is at the expense of frequency resolution, or the ability to discern between two closely spaced tones that can be approximately the same power level (with decreased resolution, they would appear as one tone).

Note that the window will reduce the total power of the signal. For the measurement of power in bins where there are signals of similar power in adjacent bins (such as noise or waveforms that do not occupy single frequency bins), the non-coherent gain of the window would be used to compensate the measurement). For measurement of a single tone that occupies one bin (such as done with a calibration measurement), the coherent gain of the window would be used. The difference (in dB) between the coherent gain and non-coherent gain is $10Log10(RBW)$, where RBW is the resolution bandwidth. This makes sense as with a rectangular window (no further windowing other than selecting the data) the RBW is 1 bit, but any other window increases this bandwidth, thus if you sum square the bins to get the total power in the waveform, if you windowed other than the rectangular window you would be over-counting if there is similar power in adjacent bins since each bin includes the power in adjacent bins. The second link below goes into that further.

More on windowing:

How to apply Hamming Window?

How to calculate resolution of DFT with Hamming/Hann window?


Product of $L^2$-norms Let $H(e^{j\omega})$ represent the frequency domain response of a minimum phase filter such that $H(e^{j\omega})^{-1}$ exists. Is there a closed-form expression or a bound for the following expression?
$$
\left\lVert H \right\rVert_2^2 \left\lVert \frac{1}{H} \right\rVert_2 ^2
$$
EDIT:
$L^2$-norm is defined as a summation over all frequencies. Obviously, the lower bound will be 1 (if averaging), and equality holds if and only if $H(e^{j\omega})=\text{constant}$. But is there a better bound for a general filter response?
Thanks.

But is there a better bound for a general filter response?

No. The argument about it $H$ being a minimum phase filter doesn't make a difference. In essence you are asking a math question on some bound for
$$L = \sum X^2 \cdot \sum \frac{1}{X^2}$$
There is none. Let's assume a very simple signal $x = [a,1]$ . The product of the norms is
$$L = (a^2+1)(1/a^2+1) = a^2 + 2 + 1/a^2.$$
It's clear that $L >=4$ and the equation is easily solvable for $a$ for any given value of $L$. In other words: by choosing $a$ you can make $L$ as large as you want it to be, so there is upper bound.


STBC Rank Criteria Explanation I am learning about space-time block coding and I am trying to understand why the rank criteria is a good measure of diversity gain. I know that we want to maximize the distance between codewords $X_i, X_j$, and I understand that $G_{ij} = (X_i - X_j)^H (X_i - X_j)$ is "bigger" when the two codewords are further apart, but i don't see how $rank(G_{ij})$ measures this matrix's "bigness."

Why does rank measure codeword detection error? It seems like we would want to consider the size of the values in $G_{ij}$, but this measure doesn't do that. Is there some eigenvalue stuff going on behind the scenes that I'm not seeing?

Thanks for any and all help!

[First shot] What I like with SE.DSP is that you can try to answer questions on topics you know nothing about. And I know nothing about space-time block coding. But I'll dig into the papers to be more precise later, since I am quite interested in concepts like diversity, redundancy, sparsity.

Most signal processing algorithms reside in the gray area between: "what you really want" and "what you can reasonably achieve". Sometimes, what you want is not computable, or difficult to model, and most of your DSP  work consists in finding a "good-enough" proxy to satisfy both goals. This the modelling and optimization part of DSP.

Concepts like diversity, redundancy, sparsity are rather vague, and could formalized in many and very different ways, depending in the objective. One acceptation of diversity is: how good am I at describing the richness of the space of possibilities? Good may means: do I describe the space of possibilities as much as I can (with respect to my dimension) or do I describe this space efficiently?

Let us be more concrete. In a 2D linear space, I only need two basis vectors to describe all other vectors, by linear combination. Wait, not any tow vectors. Two that are not co-linear. Two (and more vectors) that are colinear, or linearly dependent, have lower rank.
And we do not care about the size of the vectors. So basically, the rank is one mathematical object that both:

check how much a set of vectors is not colinear
does not really bother about their actual size (if $\vec{a}$ and $\vec{b}$ are not colinear, $\vec{a}$ and $10^{-6} \vec{b}$ ae not as well)

So, the issue is not about being big (as scalar factors barely matter), but being "spread" (or dense).

I'll be back latter, building on parts of: Space–Time Codes for High Data Rate Wireless Communication: Performance Criterion and Code Construction,  IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 44, NO. 2, MARCH 1998


Resources for an introduction to computer vision Does anyone have any suggestions on good resources or sites for learning about computer vision?

If (from the comments) you are interested in an introduction to computer vision (for the difference, see Question on difference between CV and IP), you could take a look at the introductory course of Coursera on Computer Vision..

I am also currently enrolled, but have not had much time to look into it, but it seems to me as a decent course (with the knowledge of computer vision I already have) to get an overview of the different subjects in computer vision, and the techniques that are used to overcome problems.

Next to this, there is a free download of the book Computer Vision by Richard Szeliski, which is also a great source for the basic principles of computer vision, accompanied with a lot of examples.

With this basic understanding and feeling with the subjects of computer vision, you will most likely understand more complicated things easier, because you can place things in perspective.


Discrete version of pillbox function This is an image processing question.

In 2-D Euclidean space a pillbox function, $f(x)$ is given by:

$$f(x,y) = 1\text{  if } x^2 + y^2 < R, 0 \text{ otherwise } $$

$R$ is come constant.

If I try to approximate the top hat when I'm looking at an image matrix of points, a naive algorithm would make some of the image pixels 1 inside a given radius and the others 0 outside that radius.

However, if I'm looking at a 4 x 4 grid of points, the shape I'd approximate would look like a square if I tried to naively create a discrete top-hat.

How would I create a discrete top-hat function that accurately approximates it's continuous counterpart, even when my image matrix is small?

EDIT:  I believe the comment about the edge image pixels having a value equal to the area of the pixel contained within the pillbox is correct. However, implementing this algorithm seems difficult. A reference would be appreciated.

EDIT 2: Matlab has an algorithm for computing a discrete pillbox function in it's fspecial function (disk option). This is the algorithm I'm after, though I still need a reference or an explanation of the geometry, as the function is quite complicated.

I think you should first construct a smooth function of $r$ (where $r=\sqrt{x^2+y^2}$). For example, your original solution was based on the following function of one variable:

$$f_1(r) = \begin{cases} 1 & r < R \\
0 & \text{otherwise} \end{cases}$$

but you might try using a "softer" step function such as:

$$f_2(r) = \begin{cases} 1 & r<R-\frac{s}{2} \\
0 & r\ge R+\frac{s}{2} \\
\frac{1}{2} - \frac{r-R}{s} & otherwise
 \end{cases}$$

and you can create smoother variants, such as $$f_3(r) = 3f_2(r)^2-2f_2(r)^3$$ since $f_3$ is simply formed by mapping the linear ramp function of $f_2$ to smoother cubic shape.

The transition width of the ramp, $s$, should be set to $\frac{1}{2}$ or $\frac{1}{\sqrt2}$ if the linear ramp, $f_2$, is used, and you may need a larger value of $s$ (say, $s=1$) if the cubic function, $f_3$ is used.

Note that, strictly speaking, what you are trying to do is sample the original pillar-box function at discrete sample-points, so if you think about the function prior to discrete sampling (let's call this the pre-sampling function, $f(x,y)$ being defined over the Cartesian plane, for all real $x$,$y$)  the problem with your original pre-sampling function was that it contained a step-discontinuity (a $0^{th}$order discontinuity), which contains large amounts of high-frequency energy (in 2D spatial frequency), so you need to remove some of this high-frequency energy.

The radial function, $f_2(r)$, removes the $0^{th}$order discontinuity, resulting in far less high-(spatial)-frequency energy in the pre-sampling signal (but the resulting signal still has a $1^{st}$order disconinuity, since it's first derivative is not continuous). Likewise, the cubic function, $f_3(r)$, has even more reduced high-frequency energy, since it has no $0^{th}$ or $1^{st}$order disconinuities (but it's second derivative is discontinuous).

In general, if you make sure your pre-sampling function has it's first $n-1$ derivatives being continuous (and it's $n^{th}$ is disconinuous), then the amplitude of it's high frequencies will fall of like $\left(\frac{1}{f}\right)^{n+1}$

I hope this helps...


Estimating a Signal Given a Noisy Measurement of the Signal and Its Derivative (Denoising) I have a signal and its derivative simultaneously measured, both including additive noise. The measurement is completed before the analysis, so it can be looked ahead. Now I want to reconstruct a less noisy version of the signal. I'm looking for pointers to algorithms I should look into.

Kalman filter seems to be on the right track, but the implementations I see so far are trying to estimate based on previous measurements only, while I probably should use both previous and coming measurements for optimal results at each point.

Ideas?

This is a really nice problem.
Problem Formulation
I will formulate it as following:

Let $ x \in \mathbb{R}^{n} $ be a signal. Given $ y \in \mathbb{R}^{n} $ which is a noisy measurement of $ x $ such that $ y = x + v $ and $ z $ be a noisy measurement of the derivative of $ x $ such that $ z = F x + w $ where $ F $ is the finite differences operator. It is known that both $ v $ and $ w $ are White Noise which is independent form each other. Find a reasonable estimator of $ x $.

OK, actually, I did the hard part by formulating the problem the way I did (Formulating a connection between $ x $ and its derivative).
I find this formulation to be reasonable. It doesn't assume too much but gives us enough to work with.
So my first of choice would be something like:
$$ \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| F x - z \right\|}_{2}^{2} $$
This is really easy to solve.
The tricky part is $ \lambda $. Well, I put it there to allow playing with the ratio between the noise level of the measurements. If the $ \sigma $ of both noises is the same, set $ \lambda = 1 $. Else, set it as $ \lambda = \frac{ {\sigma}_{v}^{2} }{ {\sigma}_{w}^{2} } $. The intuition is straight forward, if the noise level of $ w $ is smaller then $ v $ we want $ \lambda $ to be bigger than $ 1 $ in order to maximize the knowledge coming from $ z $.
Remark
The model above is actually the Maximum Likelihood Estimator of the model above if you assume the noises are AWGN.
I don't want to get into the Math of deriving it, but Least Squares is almost (Always) the ML for the AWGN model.
Solution to the Model
As mentioned above:
$$ \hat{x} = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| F x - z \right\|}_{2}^{2} $$
This is a strictly convex problem hence the solution is given in the stationary point.
$$\begin{align*}
\frac{\mathrm{d}}{\mathrm{d} \hat{x}} \frac{1}{2} {\left\| \hat{x} - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| F \hat{x} - z \right\|}_{2}^{2} & = 0 && \text{Definition of Stationary Point} \\
& = \hat{x} - y + \lambda {F}^{T} \left( F x - z \right) && \text{} \\
& \Leftrightarrow \left( I + \lambda {F}^{T} F \right) = y + \lambda {F}^{T} z && \text{} \\
& \Leftrightarrow \hat{x} = {\left( I + \lambda {F}^{T} F \right)}^{-1} \left( y + \lambda {F}^{T} z \right)
\end{align*}$$
MATLAB Simulation
The model I chose is using an harmonic signal - sine.
We use Finite Differences as the Derivative Model.
Here are the signals and the estimation using $ \lambda = \frac{ {\sigma}_{v}^{2} }{ {\sigma}_{w}^{2} } $:

It can be seen that the estimation is much better then any of the given signals.
Let's verify the optimality of $ \lambda $:

It is hard to see 2 circles at the bottom but they are close to verify the optimality of $ \lambda $.
The code is available on my StackExchange Signal Processing Q52150 GitHub Repository.


How are these 2 frequencies aliases of each other? I'm trying to understand this solved example from a book, given an analog signal $x(t) = 3\cos(100\pi t) $. This signal has a frequency of $50hz$. If sampled at $F_S=75Hz$ ,I obtain the discrete time signal $3\cos(\frac{2\pi n}{3}) $.

Now I am supposed to find the frequency $F $ lying in the range $0 < F < F_S/2 $ that yields identical samples. Now I get that this frequency will be $F = 25Hz$ by using $f = \frac{F}{F_S}$.

However, I learnt earlier that this frequency is supposed to satisfy the relation $F_k = F_o +kF_S$ for integer values for k, and as it is obvious, there seem to be no integer solutions for this. Where am I going wrong ?

Any help would be appreciated.

Book Reference : https://engineering.purdue.edu/~ee538/DSP_Text_3rdEdition.pdf

Note that the cosine has a positive and a negative frequency component. Consequently, in the range $[0,f_s/2]$ you get a component at $-50\textrm{ Hz}+75\textrm{ Hz}=25\textrm{ Hz}$.


Edge Detection on a Color Image I understand the process of using a Sobel kernel for edge detection in greyscale images. The input is a greyscale image, and the output is a greyscale image. I'm having trouble, however, figuring out how to apply the Sobel kernel to a color image.

Do I find the color gradient for each color channel (R, G, and B) then average each value to get the final color gradient value? Or do I average the color channels then find the color gradient of that image?

My goal is to render a raytraced image, find the edges, then re-render only the high-contrast regions with anti-aliasing and leave the low-contrast regions alone.

Finding edges in a color image can be done by decomposing the image into its channels, finding the gradients separately and fusing them somehow. However, such approach doesn't incorporate the color components in a joint model. Luckily, there is a better way to do this, which is the structure tensor representation.

The color structure tensor describes the bi-dimensional ﬁrst order differential structure at a certain point in the image. It is specified by:

$
   S=
  \left[ {\begin{array}{cc}
   R_x^2+G_x^2+B_x^2 & R_xR_y+G_xG_y+B_xB_y \\
   R_xR_y+G_xG_y+B_xB_y & R_y^2+G_y^2+B_y^2 \\
  \end{array} } \right]
$

where subscripts denote the partial derivatives. This is a more precise description of the local gradients. Eigen-decomposition is then applied to the structure tensor matrix $S$ to form the eigenvalues and eigenvectors. The larger eigenvalue shows the strength of the local image edges (gradient magnitude) and the corresponding eigenvector points across the edge (in gradient direction). In other words, eigenvalues encode the gradient magnitude while the eigenvectors contain the gradient orientation information. Note that the eigenvectors and eigenvalues can be computed analytically for a structure tensor.

There are many other usecases of the structure tensor, so it is good to know. Last but not least, here is a MATLAB code, which uses this information to compute the gradients:

http://www.mathworks.com/matlabcentral/fileexchange/28114-fast-edges-of-a-color-image-actual-color-not-converting-to-grayscale/content/coloredges.m

After obtaining the gradients, you could further carry on your edge detection in a standard fashion (canny etc.).


Amplitude of frequency in MATLAB FFT I am trying to extract amplitude of specific frequency in Matlab FFT. Is it possible to use the abs(mag)... But I do not know in which sample to look for mag(245) should give me amplitude for the frequency of that sample...

How to extract that amplitude-magnitude for $120$ Hz using mag?

I will add simple code:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise
f=(0:L-1)*Fs/L;
x=fft(y);
mag=abs(x);
mag(1)=0;
plot(f(1:L/2),(2/L*mag(1:L/2)));
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

The formula you are looking for is

  (desired freq / Sampling Frequency) * Length of samples = sample number

you can see this works out even by verifying the units

  (Hz / Hz )* sample = sample

depending on your exact frequencies, the result may or may not be an integer number. If it is not an integer, that means the exact frequency wasn't captured (due to Fs). For example if your result was 23.8 you would have to look at either sample 23 or sample 24 to get the data you required, or do some sort of interpolation.

Good luck

edit

as your comment specified, this doesn't work exactly, the reason is becasue matlab does 1 indexing, not 0 indexing. So our DC component 0Hz is at element 1, not element 0 but the above formulas are correct. However implementing in matlab is a little different

desired_freq = 120; %Hz
desired_index = (desired_freq / Fs) * L + 1; %+1 is because matlab uses 1 indexing not 0 indexing

again desired index may not be an integer number so you can't simply say mag(desired_index) one workaround would be mag(round(desired_index))


Negative group delay and envelope advance I am having a doubt reading about delays in signal processing.

Let there be an input to a LTI system with frequency response $H(f)$, given signal $x(t) = a(t)\cos(2\pi f_ot)$, where $a(t)$ is a narrowband envelope and of bandwidth $B$ Hz, such that $f_o \gg B$.

The output $y(t)$ is given by $$y(t) =\lvert H(f_o) \rvert a\left(t - \tau_g\right)\cos\left(2\pi f_o\left(t - \tau_p\right)\right)$$ here $\tau_g$ is the group delay and $\tau_p$ is the phase delay at frequency $f_o$.

My question

It is not unusual to have positive gradients in the phase response and hence have a negative group delay which seems to suggest that the envelope is advanced in time, so the input appears at the output prior to being applied!! This of course should not happen in practice so what are we missing here. Can somebody explain this. Is there an issue with the derivation? But this is a well known equation.

Note:I have read the research article and other questions  that conclude with the notion of "filter being able to predict values from past" , I am not convinced by those. In a causal LTI practical system I am sure, an input will appear at output only after it is applied at input.

Answer : No, any causal LTI system with frequency response $H(f)$ cannot produce the output $y(t)$ in advance. And, the answer lies in the causality of input signal $x(t)$ being applied to $h(t)$. Any causal input $x(t)$ which has an identifiable beginning cannot truly be Narrow-Band or Band-Limited. It will have non-zero frequency content at all frequencies.

Yes, you are right that it is not uncommon to have practical LTI systems with $+ve$ gradient of $\angle H(f)$ in parts of the response and hence making group delay $-ve$ around those parts of $H(f)$. And, if we can give a Narrow-Band input $x(t)$ such that the bandwidth of $x(t)$ is restricted in that part of $H(f)$, then you would have a time-advanced output. So, are we able to look into the future?

NO!!!! We are not. My point will get clear in a minute.

Let me take an example of a very common and practically realizable IIR filter in equivalent discrete time scenario : the Leaky Integrator.

The $H(e^{j\omega})$ of a leaky integrator is given by the following:
$$H(e^{j\omega}) = \frac{1-\lambda}{1-\lambda e^{-j\omega}},$$So,
$$|H(e^{j\omega})|^2 = \frac{(1-\lambda)^2}{1 + \lambda^2 -2\lambda cos(\omega)}, \angle{H(e^{j\omega})} = arctan \{\frac{-\lambda sin(\omega)}{1-\lambda cos(\omega)} \}$$

The shape can be plotted in MATLAB by following:

  freqz(0.1, [1 -0.9], (-pi:0.001:pi));

Now, if we can give a very narrow-band input $x[n]$ centered around $\omega = 0.6\pi$ and bandlimited within a very small $\Delta \omega$, then we would get a response as follows:

$x[n] = s[n]cos[\omega_o n]$, where s[n] is a narrowband baseband signal and $\omega_o = 0.6\pi$ and group delay of the filter is $g_d$ around $\omega_o$
$$Y(e^{j\omega}) = X(e^{j\omega}).e^{-j.g_d(\omega-\omega_o)},$$ You can work this out to get $y[n] = s[n - g_d]cos[\omega_o n]$

According to the equation above, leaky integrator is basically producing an output which is having a delayed envelope of input by $g_d$ samples. And, what happens if this $g_d$ is negative!

Check out that $g_d$ is indeed negative around $\omega_o = 0.6\pi$. Does that mean that the leaky integrator is able to produce the $s[n]$ envelope $g_d$ samples in advance?

No, it is not. The caveat is that we cannot have a perfectly bandlimited narrowband causal input $x[n]$. We cannot have a $x[n]$ which has an absolute start in time and yet it is having a bandlimited narrowband response in frequency domain.

Because we cannot have such input $x[n]$, hence we cannot have a "future seeing time machine".

In order to produce a causal input, which has an identifiable absolute start in time, the frequency response of the input will spread in frequency domain and the input $X(e^{j\omega})$ will be present at all frequencies with non-zero spectral components, and this will make the overall delay to be positive.

Indeed, if you plot the group delay response of the leaky integrator, you get the following, and check that even though the group delay is small negative number away from $\omega = 0$, it is taking high $+ve$ values around $\omega = 0$:

Hope that answers your question.


4-phase FIR LPF in MATLAB I'm trying to create a low pass FIR filter using the coefficients I found in ITU 1770-3 on page 18. This is the last piece of the puzzle for an audio loudness metering algorithm I'm working on.

I'm not super familiar with digital filters in general, but so far I've been able to create them using numerator and denominator coefficients. From there I can use filter(b,a,input) to process the audio.

But I'm not quite sure what to do with these 4 sets of coefficients. Am I supposed to get a set of numerators/denominators out of them somehow? What can I do to create an LPF out of these coefficients?

UPDATE: With some sage advice from @DanBoschen I have modeled the filter. Here is the response chart:

These appear to be the coefficients of a linear phase (due to the even symmetry of the coefficients) 48 tap low pass filter arranged in polyphase structure and given the words interpolator in the title, for purpose of upsampling the input signal by a factor of 4.

I have detailed the approach to designing a polyphase interpolator from the FIR coefficients in this post: How to implement Polyphase filter? ,  with the image copied below that specifically shows how a traditional FIR structure is mapped to a polyphase structure for interpolation using to row to column mappings of the coefficients:

In the case of this specific 48 tap implementation in this post given as four 12 tap polyphase filters, the signal would first be upsampled with a zero insert just as done in the diagram I posted above (so insert 3 zero's in between each input sample resulting in a 4x increase in the sample rate).  The low pass filter is designed to pass the signal of interest and reject the images that the zero-insert creates. (I have more details on that specifically but that's more detail than necessary to answer the question so will lean slightly toward brevity--- if anyone wants that added, please indicate in the comments). In this case given the completed design, a 48 tap filter was determined to be sufficient by the designer, and this was then mapped into four 12 tap poly-phase filter structures (instead of four 2 tap filters in my image above) as follows:

For FIR filter with taps h[0], h[1], h[2].... h[47]

This is mapped into four poly-phase filters as:

Filter 1:
h[0] h[4] h[8] h[12] ... h[44]

Filter 2:
h[1] h[5] h[9] h[13] .... h[45]

Filter 3:
h[2] h[6] h[10] h[14] .... h[46]

Filter 4:
h[3] h[7] h[11] h[15] .... h[47]

The columns given by the OP are the taps for each of the four filters as detailed above. To view the filter response, read the coefficients back in order from each row as listed and then can use the freqz command available in Matlab, Octave or Python (SciPy) to view the response.


Can anyone suggest good algorithms for CBIR? Project: Content Based Image Retrieval - Semi-supervised (manual tagging is done on images while training)

Description

I have 1000000 images in the database. The training is manual (supervised) - title and tags are provided for each image. Example: coke.jpg Title : Coke Tags : Coke, Can

Using the images and tags, I have to train the system. After training, when I give a new image (already in database/ completely new) the system should output the possible tags the image may belong to and display few images belonging to each tag. The system may also say no match found.

Questions:

1) What is mean by image fingerprint? What is the image fingerprint size expected ? (important because there will be millions of images to be inserted in database)

2) What is the field format of that fingerprint in the database ? (important because a fast search is needed … script should search in a 1M images database in less than 1 second)

3) What is the descriptors (algorithms) we use to analyze them ?

Thanks in advance

A luminosity histogram (especially one that is separated into RGB components) is a reasonable fingerprint for an image - and can be implemented quite efficiently.

Subtracting one histogram from another will produce a new historgram which you can process to decide how similar two images are.
Histograms, because the only evaluate the distribution and occurrence of luminosity/color information handle affine transformations quite well.
Luminosity histograms produce false negatives when the color information in an image is manipulated.


Problem with calculating recursive equation I'm facing with this error when I'm trying to solve the equation :

  Attempted to access x(101); index out of bounds because
  numel(x)=100.

actually the original equation is

so I used "i+1"instead of"i" to avoid following error:

  Attempted to access y(0); index must be a positive
  integer or logical.

n=[1:100];
yy=zeros(1,99);
xx=zeros(1,99);
y=[1 yy];
x=(.5).^n;
for i=n

  y(i)=.5*y(i-1)+2*x(i)+4*x(i-1)
end

error resolved by initializing i by  2

the improved code is :

yy=zeros(1,99);
xx=zeros(1,99);
y=[1 yy];
x=(.5).^n;
 for n=2:100

 y(n)=.5*y(n-1)+2*x(n)+4*x(n-1);
end

but the results of **y(n)**is different from what I get on paper.

Sorry for answering like this but I cannot comment since my reputation is not sufficient, but yet I want to help.

It is a filter and you can do it way easier, and more quickly than just building an algorithm(in Matlab ofcourse). It will take relatively a long post to cover all the concepts in order to make a sufficient understanding about filter designs.

However, I can help you about which topics you should cover.

Please take a look at FIR/IIR filters (Your recursive equation is basically an IIR filter) and Z-Transform topics.

If you need a quicker way to learn about IIR filter this lecture will give you some information.


Second order sections and normalized factor I use matlab function tf2sos to convert digital filter transfer function data to second-order sections form:

  [sos,g] = tf2sos(b,a)

Taken from the help menu of the matlab:

  G is a scalar which accounts for the overall gain of the system

My question is: Because G is the overall gain if i draw the block diagram is going to be something like below (i use office visio for the pictures, if you press view picture the resolution is good but in this thumbnails they look awfull, sorry about that) :

But if i want yo design my system like:

how can i find the gain between every node(b1,b2...etc)?

If I understand your question correctly, you want to split the total gain between two or more stages. In this case you just have to satisfy $G=g_1g_2g_3$, where $G$ is the total gain returned by tf2sos() and the $g_i$ are the individual gains per stage.


How is RMS and Peak-to-peak of white noise calculated I have seen some devices which specify their internal noise in RMS or Peak-to-peak value. How does one calculate RMS or p2p values for white noise?

three things to think about:

"white" describes the power spectrum, but to know the peak values, you need to know the probability density function (p.d.f.) of the noise process.  often, with "white noise", we are assuming "normal" or "gaussian" p.d.f.
even for a finite power (which means a finite variance for the p.d.f.), a gaussian random variable has infinite peak value (just that it hits whatever arbitrary high value with decreasing probability).  so for gaussian, you need to define like a "1% likelihood peak value".
true white noise (which does not exist in reality) has infinite power and infinite variance.  this is because true white noise has infinite bandwidth.  so, in reality, we deal with bandlimited "white" noise, so the bandwidth might be a necessary parameter to know, before you can think about its r.m.s. (which is the square root of the variance, assuming no DC component) and peak.


how to extract meaning from signal (distance estimation)? I have 3 times series representing detected rssi signal power from 3 emitting devices. The devices are at 1 meter distance from the receptor, plotting the 3 time-series gives the following results :

My question is : Is there a way, a general mechanism, a filter, or an algorithm for processing these signal, so that I can find some common regularities between them to which I can assign the distance 1 meter, and predict that distance for different time-series?

In general, you can apply a filter with some time constant T which will be giving you a more reliable RSSI value but even in this case, the RSSI will only be an accurate indication of distance under certain assumptions.

There are two types of filters you could use here, a moving average low pass filter which would effectively produce a smoothed out version of the RSSI waveform or a non-linear modal-mean or median filter. An example of an easily implemented moving average Finite Impulse Response (FIR) filter would be h=ones(N)./N where N is the order of the filter that is directly proportional to the Group Delay (or time constant) of the filter (for instance a [0.25,0.25,0.25,0.25] N=4 filter). Applied over a sliding window of N samples, this filter would return an average, slow changing, estimate of the RSSI curve which you could then associate with distance. Similarly, the non-linear filter is applied over a sliding window of size N but in this case you just estimate the modal mean (i.e. the most frequently occurring RSSI value of the N available in the sample) or the median value. The median filter would tend to remove the spikes and could possibly be a better alternative for this application.

The relationship between RSSI and distance would be established through the Free Space Path Loss but the main assumptions here would be that a) there is a direct line-of-sight between transmitter and receiver and b) the RSSI is accurately estimated by the device. In other words, if you start moving the device in space, it is not necessary that the distance estimation will remain accurate. You might also want to check this link as well.

Hope this helps.


Calculating the cutoff frequency of a highpass digital filter I am new to DSP and I am trying to find the cutoff frequency of a HP digital filter. I know the equation that describes the system, its frequency response: $H(e^{jω})= 1 - \frac{e^{jω}}{2} - \frac{1}{2e^{jω}}$ and of course amplitude response / phase diagrams.

I have found (p. 8) that for analog filters I can calculate $\frac{1}{\sqrt{2}} |H(e^{jω})|$. However I don't know how to do this for digital filters. Reading this makes it look like it's the same procedure, although this seems (to me) like a contradiction to

  Digital filters are less standardized, and it is common to see 99%, 90%, 70.7%, and 50% amplitude levels defined to be the cutoff frequency.

of the first source. Am I confusing something? Please keep in mind that I am new to DSP. Thank you.

Concerning the cut-off frequency, there's really not much of a difference between analog and digital filters. $3\textrm{ dB}$ is common, but any other value is fine, as long as you specify it and people know what you're talking about.

In the case of the given filter, it's quite straightforward to compute any cut-off frequency. Let's choose the $3\textrm{ dB}$ cut-off frequency:

$$H(e^{j\omega})=1-\frac12\left(e^{j\omega}+e^{-j\omega}\right)=1-\cos(\omega)\tag{1}$$

Since $H(e^{j\omega})$ is real-valued and since $H(e^{j\omega})\ge 0$, we can simply solve the equation

$$1-\cos(\omega_c)=\frac{2}{\sqrt{2}}=\sqrt{2}\tag{2}$$

Note that the maximum of $H(e^{j\omega})$ is $2$, so the value of $H(e^{j\omega})$ at the $3\textrm{ dB}$ cut-off frequency $\omega_c$ is $\sqrt{2}$.

From $(2)$ we get $\omega_c\approx 1.9979$. The actual cut-off frequency in Hertz is given by

$$f_c=\frac{\omega_c}{2\pi}f_s\tag{3}$$

where $f_s$ is the sampling frequency.


Estimating a phase for a frequency band I'm looking for a way (preferably already implemented in eg. python) of estimating a phase for a specified freq band.
When I'm performing FFT on the signal, I can determine components frequencies as well as phases of interesting components. If I understand this well, a single complex number in the FFT result table is, in fact, a representation of some (small and determined by the sampling rate and length of signal) freq band. But what if I want specific a frequency band on my own? Eg. I want to estimate a phase (single number) for a frequency band of 9-11 Hz. I guess even it's possible, some information will be lost. But can I do that?
Thank you in advance.

a single complex number in the FFT result table is, in fact, a representation of some (small and determined by the sampling rate and length of signal) freq band

Not really. The DFT is first and foremost a mathematical procedure that simply allows representing a time domain sequence as a sum of complex exponentials. The physical interpretation of that operation depends a lot on the circumstances but in general it's NOT true that FFT coefficient at some frequency $\omega$ represents accurately what's happening physically at frequency $\omega$. They are related but the exact nature of the relationship depends a lot on the details.
The notable exception here are time sequences that are periodic with the FFT length, but periodicity implies infinite length which doesn't exist in the real world.

want to estimate a phase (single number) for a frequency band of 9-11 Hz.

That would require you to first define what a single number phase for a frequency band is. Typically the phase is different at each frequency, so typically a single number doesn't
work well.


Determing inverse Z-transform using impulse response? In Matlab there is a command iztrans for finding inverse Z-transform. But how can we find inverse Z-transform using impulse response? The Matlab command impz gives Impulse response of digital filter but how can I then proceed towards finding inverse Z-transform?  Actually i want to determine inverse ztransform of H(z) shown in attached photo. But i want to determine inverse z transform using impulse response since i read somewhere that inverse ztransform can be determined by impulse response using matlab but i am unable to find any google resource that guides much in this regard. I have also attached a table of common z tranforms. I have also highlighted the two most relevant cases as far answer of GKH is concerned

When you have an $H(z)$ (a transfer function), its inverse Z-transform is the impulse response $h[n]$. iztrans is a symbolic function to compute any (causal) inverse Z-transform, either it is related to an impulse response or not. impz is a non-symbolic function to compute some samples of the impulse response of a digital filter.

If I were you and I wanted to find the inverse Z-transform of the given $H(z)$ - which is the impulse response -  I would use residuez:

[r,p,k] = residuez([4 3 9], [4 3 -4])
r =

 1.0544
 2.1956

p =

-1.4430
 0.6930

k =

 -2.2500

That allows you to write the impulse response as $$h[n] \approx -2.25\delta[n] + 1.0544(-1.443)^nu[n] + 2.1956(0.693)^nu[n] $$

given that we're talking about a causal system.

Alternatively, you can use

[h,n] = impz([4 3 9], [4 3 -4]);

but that will give you only a couple of samples of the impulse response.


How can I use time-series transforms for detecting malware? I am student of network security, and very unfamiliar with the domain of signal processing. Do pardon my naivety at asking very basic stuff.

My work deals with detecting malicious activity in network logs -say, differentiating between botnet traffic and benign Internet traffic. I came across this interesting article titled Detecting and identifying malware: a new signal processing goal. Since then, I have been exploring the world of FFTs, DTFTs, Wavelet transforms, cepstrum and the ilk. But I have not been able to make much progress with the application of these principles.

I wish to know from dsp experts that what kind of time-series transforms can be useful for my task.

To keep things clear, I will provide a brief description of data I deal with: As I said, my task is of separating botnet traffic from Internet traffic. Lets define Internet traffic as the traffic seen from the systems of regular users of the Internet like you and me, who browse the web, watch a few videos, check mail, maybe download a few movies, and also exhibit periods of inactivity say while they sleep, eat, etc.

And lets define botnet traffic as the traffic of a covert malware/virus (a 'bot') on your machine calling home to its master or contacting fellow 'bots' (just to maintain connectivity), and also the commands sent by the master to its bots or those exchanged by bots amongst themselves. This traffic does typically have patterns and periodicity.

If I am right, I am dealing with signals of botnet traffic mixed with Internet traffic (which does not have any specific pattern?) What kind of time-series transforms can be useful here? DFTs? FFTs? Cepstral analysis {or is that being too creative? :)}

The point seems to be moot. An attacker could hide in noise. Any meaningful interpretation of the question appears to aim for temporal regularity of attacks. Maybe an attack is temporally regular, but so is users' behavior. So, there might be lots of applications for transforms in malware detection, but none of them appears to promise any insight on how to recognize or intercept it.

From a computer science standpoint i would say that ,,recognizing malware'' is akin to the problem of recognizing two functions as equal, which in general is a proof, which is not computable in any ordinary sense. It would probably amount to proving P=NP. Any ,,heuristic'' recognition of malware can be evaded that way. The only system capable of actually and reliably recognizing malware appears to be a conscious human being.


Frequency of a single sample in a digital signal and aliasing I'm involved with digital audio synthesis. I know that if I create a raw non band-limited waveform it would contain frequencies above Nyquist thus violating the Nyquist Shannon theorem and the signal will contain aliased components. So I wonder if I can also create aliasing by other means, like a "click" from a fast envelope or a single sample in a digital stream. So what is the frequency of a single sample (is it sample rate/1 ?) and does such clicks alias in the frequency domain?

A single sample being non-zero is known as an impulse.  The effect in the discrete frequency domain depends upon where it is within a the frame.  The DFT will give you the coefficients required to reproduce whatever finite signal is withing the sample frame.  It is misleading, even erroneous, to say the signal is composed of those frequencies.  This can only be said about periodic signal for which the sampling frame is a whole number of periods.  Recall that the DFT and inverse DFT are basically the same, except for a sign change in the exponent.  When you have a pure tone with a whole integer frequency in the sample frame, the result is a single spike in the DFT at the bin corresponding to that frequency.  So when you have a single spike, or impulse, in the signal domain, the result will be a DFT with the bin values being the same as a pure complex tone.

"Alias" is a different concept entirely.  It has to do with the sampling of tones higher than the Nyquist frequency appearing to be tones of a lower frequency.  The easiest way to envision what alias frequency will look like is to imagine the bins of the DFT arranged along the unit circle.  The Nyquist bin will occur at $\pi$ radians.  Any alias frequency has to appear somewhere on the circle.  It may be beyond $\pi$, or even beyond multiple times around the circle.  The frequency it will appear to be (for real valued signals) will be the frequency associated with the bin number of the first half of the circle that has the same projection onto the horizontal axis.

I would recommend that you read my blog articles (link on my profile page) from the start to get a perspective on this material.

Hope this helps.

Ced


What are the rough steps in obtaining and processing an image from an infrared camera? I have asked the same question in stackoverflow, but since I am not getting any responses , I am going to ask it again here.

My project is to design an interface between an IR camera (Xenics XTM-640) and an ARM/FPGA unit and then process the data to make some comprehensible image.

At the moment I am working on creating timing waveforms to control sensor and capture of frames and I am planning on how to go on from there onwards.

From what I can see there is number operations which can be performed to make the image clearer. Terms I have came about such as:

quantisation
convolution
Gaussian filterring
contrast enhancment
noise cleaning/noise smoothing
neighbourhood operation
normalisation
edge crispening

What is the difference between preprocessing and processing?
Is there a preferred order to performing this steps?
Are some absolutely necessary, while others not so much and can be skipped?
How would this differ between processing normal light vs IR radiation?

Ideally I've been looking for a flowchart/diagram that shows all (or some) of the stages from capturing IR radiation to displaying false colour image on a screen.

In the "Digital Image Processing" by William K. Pratt it seems the chapters are broken down into

CONTINUOUS IMAGE CHARACTERIZATION
DIGITAL IMAGE CHARACTERIZATION
DISCRETE TWO-DIMENSIONAL LINEAR PROCESSING
IMAGE IMPROVEMENT
IMAGE ANALYSIS

Is that roughly the correct chronologically description of the whole process?

It's hard to give a lot of useful information without knowing the details of the signals you are getting from your hardware. Typically the output from any camera is an image, so you would need to decide what is wrong with the image before trying to come up with image processing techniques to solve that problem, you don't seem to have a clear problem specification at the moment.

  What is the difference between preprocessing and processing?

Semantics. It doesn't refer to anything specific in the signal processing domain.

  Is there a preferred order to performing this steps?

Yes but you need to decide which are necessary before worrying about the order.

  Are some absolutely necessary, while others not so much and can be skipped?

It depends entirely on your source data. If your IR camera outputs a clear crisp image that contains all of the information you need then none of the steps will be necessary. If the image is completely washed out then you will need to do contrast enhancement, otherwise it can be skipped. Really nothing useful can be said without more information about your data.

  How would this differ between processing normal light vs IR radiation?

The usual caveat applies that it depends on your data but in principle there is no difference, you are just capturing a different part of the light spectrum.

  Is that roughly the correct chronologically description of the whole process?

No not really, the process will be informed by your source data and your desired output. It's unlikely that all of those topics will apply, more likely one or two of them will and not necessarily in any particular order.

If you want more concrete advice I suggest you explain your procedure in more depth, with the specifics of the data you are acquiring and the results you expect.


Authoritative DSP book/reference material on Psychoacoustics I am an auto didactic Electrical Engineer currently pursuing a goal to design a custom application targeted at deeply embedded linux(think BBB or even Raspberrypi-2) type device. My DSP background is mostly theoretical and dates back to my college days. I am currently referring to Understanding Digital Signal Processing- Richard Lyons. A tome it definitely is, but is a very basic reference at the best. I want more concrete book/reference material that helps me understand the various algorithms, their trade-offs and implementation details(pseudo code or Matlab code would definitely help me save time). For a beginner such as me who has no industry bearings, terms such as YIN,YAPPT,HPS,Cepstrum,ASDF,AMDF and many more TLAs are mere jargons. I want to do a deep dive with whatever DSP basics i have and complement that with the said reference book.

My plan is to investigate several algorithms, simulate them in Matlab to verify results and after selecting one based on the results of my experiments move on to actual implementation. My end requirement is that the Algorithm
 has to be real-time. I basically want to hook up a Beaglebone black to my son's piano and be able to detect Pitch and onset in real-time. I have been reading some interesting questions posted here and elsewhere.But I need a strong footing on Psychoacoustics  before any of it starts making sense to me.

The Fastl-Zwicker book is usually the introduction to this topic.  The more I learn about it the more fascinating it becomes.  If you pick it up used be sure to get one with the CD examples in it, as they are completely fascinating.


change function from continuous to discrete I have function that I want to sample with sampleing rate 20 Msamples/sec.

The funciton is the one below

$w_T(t)= \left\{ \begin{array}{ll}
         \sin^2(\frac{\pi}{2}(0.5+\frac{t}{T_{TR}})) & \mbox{if $-T_{TR}/2 < t<T_{RT}/2$};\\
        1 & \mbox{if $T_{TR}/2 < t<T-T_{TR}/2$}\\
 \sin^2(\frac{\pi}{2}(0.5-\frac{t-T}{T_{TR}})) & \mbox{if $T-T_{TR}/2 < t<T+T_{TR}/2$};\\
\end{array} \right. $

I saw the solution is

$w_T(n)=\left\{ \begin{array}{ll}
         1 & \mbox{if $1<n<79$};\\
        0.5 & \mbox{if $0,80$}.\end{array} \right. $

I tried to see why it so so I put $t=n T_s= n \frac{1}{20\times10^{-6}}$. It is also given that $T_{TR}= 100 ns$ and $T=4
\mu sec$. But i fail to reach the discrete funciton.

Is there any way to reach solution?

if you take sampling frequency as 20Msamples/sec then sampling time Ts is
$$T_s=\frac{1}{f_s}\\T_s=\frac{1}{20\times10^6}\\T_s=50ns$$since the solution is till n=79 we can see that t=79*50=3950
From the continues time function this  Trt/2=50ns and T-Trt/2=3950ns hence we get 1 n=1 to 79 when n=0 we get the first case which is equal to 0.5 similarly of n=80.This can be see as$$
w_T(nT_s)=
\begin{cases}
    \sin^2(\frac{\pi}{2}(0.5 +\frac{n50ns}{100ns})),& \text{if } -50\times10^{-6}<n50\times10^{-6}<50\times10^{-6}\\
    1,              & \text{if  } 50\times10^{-6}\leq n50\times10^{-6}\leq 3950\times 10^{-6}\\
 \sin^2(\frac{\pi}{2}(0.5 -\frac{n50ns-4000ns}{100ns})),&\text{if } other case
\end{cases}$$
which comes down to
$$
    w(n)=
\begin{cases}
    \sin^2(\frac{\pi}{2}(0.5 +\frac{n}{2})),& \text{if } -1<n<1\\
    1,              & \text{if  } 1\leq n\leq79\\
 \sin^2(\frac{\pi}{2}(0.5 -\frac{n-79}{2})),&\text{if } 79<n<81  & \text{otherwise}
\end{cases}
$$
In the above equation we can remove the $$-1<n<1$$as n=0 when n=0 w(0)=0.5 for n=1to78 w(n)=1.When n=80 the third one become zero.Hence you get the finalized expression.I hope this answers your question


Homework Help: What does $h[0] = 1$ represent? What is $\ln \big| H(e^{j \omega})\big|$? I have been staring at this problem for a week now...
Suppose $H(e^{j \omega})$  is the frequency response of a stable and causal minimum-phase discrete-time system with $h[0]=1$  ($h[n]$ is the impulse response of the system). Show that:

$$ \int\limits_{-\pi}^{\pi} \ln \big| H(e^{j \omega})\big|^2 \, \mathrm{d}\omega = 0 $$

So:

Minimum phase means poles and zeros are inside the unit circle
Stable and causal is redundant since the system is minimum phase
What is conceptually meant by $h[0] = 1$?
What is the significance of $\ln$(magnitude of freq response)?

The transfer function of a discrete-time minimum phase system described by a linear difference equation with constant coefficients can be written as

$$H(z)=h[0]\cdot\frac{\prod_{k=1}^M(1-c_kz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})},\qquad |c_k|\le 1,\quad |d_k|<1\tag{1}$$

where $h[0]$ is the first sample of the corresponding impulse response, and $c_k$ and $d_k$ are the locations of the zeros and poles of $H(z)$, respectively. Note that factoring out $h[0]$ makes sure that the numerator and denominator polynomials are monic (i.e., have a leading coefficient equal to $1$).

From $(1)$, the squared magnitude of $H(z)$ is

$$|H(z)|^2=h^2[0]\cdot\frac{\prod_{k=1}^M|1-c_kz^{-1}|^2}{\prod_{k=1}^N|1-d_kz^{-1}|^2}\tag{2}$$

and its logarithm is

$$\log |H(z)|^2=\log\left(h^2[0]\right) +\sum_{k=1}^M\log|1-c_kz^{-1}|^2 - \sum_{k=1}^N\log|1-d_kz^{-1}|^2\tag{3}$$

Evaluating on the unit circle $z=e^{j\omega}$ and integrating over the interval $[-\pi,\pi]$ gives

$$\int_{-\pi}^{\pi}\log |H(e^{j\omega})|^2d\omega=2\pi\log\left(h^2[0]\right) +\\+\sum_{k=1}^M\int_{-\pi}^{\pi}\log|1-c_ke^{-j\omega}|^2d\omega - \sum_{k=1}^N\int_{-\pi}^{\pi}\log|1-d_ke^{-j\omega}|^2d\omega\tag{4}$$

The integrals on the right-hand side of $(4)$ are of the form

$$\int_{-\pi}^{\pi}\log|1-ae^{j\theta}e^{-j\omega}|^2d\omega\tag{6}$$

with real-valued $a$ and $0\le a\le 1$. With

$$|1-ae^{j\theta}e^{-j\omega}|^2=1+a^2-2a\cos(\omega-\theta)\tag{7}$$

the integral $(6)$ can be written as

$$\int_{-\pi}^{\pi}\log\left(1+a^2-2a\cos(\omega-\theta)\right)d\omega\tag{8}$$

Since the integrand is $2\pi$-periodic and we're integrating over one period, the integral is independent of the angle $\theta$, which finally gives

$$\int_{-\pi}^{\pi}\log\left(1+a^2-2a\cos(\omega)\right)d\omega\tag{9}$$

As shown in the answers to this question on math.SE, the integral $(9)$ equals zero if $0\le a\le 1$, i.e., if $H(z)$ is minimum phase. This means that if $H(z)$ is a minimum phase system, we get from $(4)$

$$\int_{-\pi}^{\pi}\log |H(e^{j\omega})|^2d\omega=2\pi\log\left(h^2[0]\right) \tag{10}$$

and with $h[0]=1$ (which is just an appropriate scaling) we finally arrive at

$$\int_{-\pi}^{\pi}\log |H(e^{j\omega})|^2d\omega=0 \tag{11}$$


Subtract audio to another audio suppose I have a sample of a snare and that snare is used in a song,¿ there is an algorithm or tool that detects the sound of the snare and subtracts it, that is to say that cancels its phase or eliminates it?.

the result would be the song without the snare sound.

this could also apply to (accordion trumpets or any other instrument)

That's a really hard problem. If you have access to the original snare stem and the amount of post processing (effects, panning, mastering) is moderate you may be able to pull some of it out with an adaptive filter.

There are a few commercial plug ins that do this type of thing, but they all have their strengths and weaknesses so it depends a bit on your specific needs


What is the "pitch marker" in the LP-PSOLA algorithm? The LP-PSOLAR algorithm uses a the Linear Prediciton Coding (LPC) to calculate an error vector from a signal. This error vector is used to calculate the pitch markers. What do these pitch markers represent in the original signal in generall? Do they represent a change in fundamental frequency or a change in phoneme (I think both assumptions are wrong)? The following figure depicts the error vector (red) and the signal (green). The second picture depicts an additional example.

The pitch marker indicates the beginning of each cycle of the waveform - or from a more physiological point of view the point at each the periodical motion of the glottis causes a sudden change in air flow.

A stable sound (in pitch and phoneme) will still have a marker every cycle - so these markers do not indicate any change.


$2\pi$ periodicity of discrete-time Fourier transform In my signals and systems course, we have learned that the discrete-time Fourier transform is $2\pi$ periodic, but the continuous-time Fourier transform is not periodic in general. For reference, we are using the following definitions of each transform:

Continuous-Time:

$$
x(t) = \int_{-\infty}^{\infty} X(j\omega) e^{j\omega t} \, d\omega \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ X(j\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt
$$

Discrete-Time:

$$
x[n] = \int_{\langle 2\pi \rangle} X(e^{j\omega}) e^{j\omega n} \, d\omega \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}
$$

I'm searching for some intuition as for why the DTFT is periodic, but the CTFT is not. In class, my instructor presented the following argument: for a discrete-time signal,

$$
e^{j\omega n} = e^{j\omega(n + 2\pi)}
$$

and thus any $x[n]$ can be expressed as a sum of individually $2\pi$ periodic functions. However, I don't see why that argument only applies to discrete-time signals - I feel as if it also works for continuous-time signals.

Any someone explain?

The argument does not work in continuous time. In discrete time the argument is that
$$e^{j\omega n}=e^{j(\omega+2\pi)n},\qquad n\in\mathbb{Z}\tag{1}$$
This is true because by definition $n$ is an integer. In continuous time we generally have
$$e^{j\omega t}\neq e^{j(\omega+2\pi)t},\qquad t\in\mathbb{R}\tag{2}$$
because $t$ is a real variable. Eq. $(2)$ only holds with equality for integer $t$, but not for all other values of $t$ (which would be necessary for $e^{j\omega t}$ to be $2\pi$-periodic).
But the fundamental reason why the DTFT is periodic is the fact that we deal with sequences, and if we assume that the sequence is the result of sampling a continuous-time signal, then it becomes obvious that different continuous-time signals can result in the same sampled sequence. E.g., take a continuous-time cosine
$$x(t)=\cos(2\pi f_0 t)\tag{3}$$
Now we sample $x(t)$ at instances $nT=n/f_s$, where $f_s$ is the sampling frequency:
$$x[n]=\cos(2\pi f_0 nT)=\cos(2\pi n f_0/f_s)\tag{4}$$
I leave it as an exercise to you to show that the frequencies $f_0=f_s-\Delta f$ and $f_0=f_s+\Delta f$ (with some $\Delta f>0$) result in the same sequence $x[n]$. This phenomenon is called aliasing. This means that frequencies above $f_s/2$ (also called the Nyquist frequency) cannot be distinguished from the corresponding frequencies below the Nyquist frequency. That's why the DTFT is periodic. There cannot be any new spectral information outside the frequency interval $[-f_s/2,f_s/2]$. If you haven't learned it yet, check out the sampling theorem, which will also help you understand the nature of spectra of discrete-time signals.
Also have a look at this answer to a related question.


How can I find expansion coefficients of the a vector in a given basis? How can I find the coefficient of the vector $\mathbf y$? And how can the inner product be done on these vectors?

Let $\mathbf y = \begin{bmatrix}1\\2\\0\\1\end{bmatrix}$
What are the expansion coefficients of $\mathbf y$ in the basis $\left\{\mathbf v^{(0)}, \mathbf v^{(1)}, \mathbf v^{(2)}, \mathbf v^{(3)}\right\}$ where
$$\mathbf v^{(0)} = \begin{bmatrix}\frac 12\\\frac 12\\\frac 12\\\frac 12\end{bmatrix}, \mathbf v^{(1)} = \begin{bmatrix}\frac 12\\\frac 12\\-\frac 12\\-\frac 12\end{bmatrix}, \mathbf v^{(2)} = \begin{bmatrix}\frac 12\\-\frac 12\\\frac 12\\-\frac 12\end{bmatrix}, \mathbf v^{(3)} = \begin{bmatrix}\frac 12\\-\frac 12\\-\frac 12\\\frac 12\end{bmatrix}, \quad $$

This is actually a math problem, but it is related to signal processing via the concept of signal space. Since it's a homework-style problem you'll only get hints that should help you solve the problem on your own.

verify that the given vectors are orthonormal.
if you've convinced yourself that the vectors are orthonormal, you can conclude that the expansion coefficients are simply the inner product of the given vector with the basis vectors, i.e., $$\mathbf{y}=\sum_{i=0}^3(\mathbf{y}\cdot \mathbf{v}^{(i)})\mathbf{v}^{(i)}$$

I don't understand your question about how to compute the inner product. If you know (or look up) the definition of the inner product then you should be able to compute it.


How to convolve $u(-t)$ with other signals? How can I convolve the following
$u(t+1)*u(-t)$
I know that convolution with $u(t)$ gives the integral of a function but what change occurs due to $u(-t)$?

In general, convolution of a function $f(t)$ with $u(-t)$ results in a function $g(t)$ defined by the following integral:

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

In the given example, $f(t)=u(t+1)$, so you get

$$g(t)=\int_{t}^{\infty}u(\tau+1)d\tau=\int_{t+1}^{\infty}u(\tau)d\tau\tag{2}$$

I leave it up to you to figure out what the result of that integral is.


Chirp z algorithm clarification I am attempting to implement a chirp z algorithm to handle random sized DFTs, and I can not seem to obtain any meaningful results. I have gone over several write ups and "think" I have a handle on what a chirp transform should do, but I wanted to ask here to double check and make sure my understanding of the process is sound. In pseudo code, here is the process I am applying.

*function Chirpz(input x[n], discrete values such that n = 0, 1, ... N - 1)*
*// To keep numbers in perspective I will assume N = 5.*

M = The minimum power of 2 which is greater than 2N - 1 // In our case of N=5 => M=16.

Pad x[n] with 0s so that it has a total length of M.

// Create y[n] by performing the following:

value = 1 - N

for n = 0 to n = 2N - 1

y[n] = e^i(PI / N * value * value)     // W^((k - n)^2 / 2)
value = value + 1

next n

Pad y[n] with 0s so that it has a total length of M.

// Multiply x[n] by W^(-n^2 / 2):

for n = 0 to n = N
x[n] = x[n] * e^-i(PI / N * n * n)

Perform FFT(x[n]) which has M elements using a radix-2 FFT algorithm
Perform FFT(y[n]) which has M elements using a radix-2 FFT algorithm

// With the results of the FFTs now in x[n] and y[n], multiply each value:
for n = 0 to n = M
x[n] = x[n] * y[n]

Perform IFFT(x[n]) which has M elements using radix-2 FFT algorithm

// Multiply x[n] by W^(-k^2 / 2)
for n = 0 to n = M
x[n] = x[n] * e^-i(PI / N * n * n)

// Now the first N elements in x[n] "should" have the result
X[k] = x[n] for n from 0 to N - 1

return X[k]

But my results don't match anything near what is returned by fftw3 given the same input. In fact my results look to be garbage. As I am pretty new to this, I just want to make sure my understanding of the algorithm is sound.

Thanks.

Ok, after some head scratching and re-re-re-reading I realize what I was missing. When creating the chirp signal used in convolution, the values of (k - n) which result in a negative need to be wrapped around and placed at the end of the signal rather than in the beginning.

So, replacing:

// Create y[n] by performing the following:

value = 1 - N

for n = 0 to n = 2N - 1

y[n] = e^i(PI / N * value * value)     // W^((k - n)^2 / 2)
value = value + 1

next n

With the following:

// Create y[n] by performing the following:

// Take care of values of k - n which are greater than zero.
for n = 0 to N - 1
y[n] = e^i(PI / N * n * n)

// Add the zero padding.
for n = N to M - N
y[n] = 0

// Now place the negative values at the end of the signal
for n = 1 - N to -1
y[M + n] = e^i(PI / N * n * n)

Yields the expected results.


Anybody with experience with online double integration of acceleration? I've found good literature for offline proccesses, and I've gotten good results. However, when it comes to online applications, I've found nothing good. Applying any filter without a linear phase response will give me too much distortion once I get to displacement. However, any filter I've used with linear phase response gives me way too much "time lag" (apart that it's very high order). Below is the proccess most people apply:

Raw Acc --> [High Pass Filter (HPF)] --> [Integrate (Int)] --> [HPF] --> [Int] --> [HPF] --> Displacement Signal

As you can imagine, this procces is trivial if you use non-causal filters... But I've not gotten a good solution for the online proccess.

My current conclusion is that I need a "low order, linear phase, causal filter." There is active research in these types of filters, but they are somewhat complicated to apply. So before I start making leeway, I decided to ask this in a forum to see if anybody could give some input or fresh ideas.

Any filter within your chain can remove useful information. The kalman filter way is applicable when you have/use a good model of the error, and you have some auxiliary sensor (like GPS, or images) to confront to the result of the double integration (DI). I dont know how good your acc is, but a military one, for example, can drift so few after the DI, that it doesnt need any auxiliary sensor. However, you cannot expect the same result using this cell-phone acc, though. The limit of goodness is defined by your application and the type of sensor you have.


What is the difference between Goertzel algorithm and cross correlation in time Im reading https://en.wikipedia.org/wiki/Goertzel_algorithm,
when we want to discover the amplitude of a some frequency in a given signal what is the difference between applying Goertzel algorithm and doing cross-correlation in time between our signal and a complex exponential with the desired frequency.
It seems like in the button line Goertzel algorithm does exactly that,
then why couldn't we apply the correlation in the first place and instead used the IIR filter ?

Thank you.

The Goertzel algorithm gets away without calculating a complex product with $e^{j\pi lk}$ n times; thus drastically reducing the computational complexity.


Sampling Theorem: How to know the value between two samples of a Signal According to Sampling theorem, in order to reconstruct a signal we need to sample it at the rate => twice the highest frequency component of that signal. (provided signal is band limited).

Let's say, we have a signal with f= 2MHz (highest freq component), so we will be sampling it at 4MHz or more as Sampling theorem say's.. and we will have N Samples in signal.

Now, what if I want to know the exact value of a sample (information) between two samples..??

How can we find the value at every possible instantaneous.. in a signal..when we know we only have N samples after sampling..??

I know we can increasing sampling rate to have more sample. But is there any other way to do it..?

This question was asked to me in an interview at NCRA-TIFR, Pune, India.

Interpolation using a Sinc function kernel.  This assumes that the sampled signal was perfectly bandlimited to below half the sampling rate.

Note that perfectly bandlimited signal are infinite in extent.  For finite-length "real world" signals, using a windowed Sinc interpolation kernel (thus a finite computation with a noise floor) is a common method to approximate points between samples (for upsampling, for instance).


How do I simulate a purely real signal being received on multiple antennae in an array? I have a signal created at baseband and then transmitted as $$\Re\{s_{bb}*e^{2\pi f_ct}\}$$ which is then received at a 5 element array for further processing. How do I interpret the real signal at each of the array elements without using complex numbers? Simulink needs a purely real signal as an input to the analog filter I am using.

The input at each array will be complex, because the wireless channel "spreads" the energy of the in-phase and quadrature parts of the signal into each other.

One simple way to think about this is transmitting a pure in-phase carrier: $$s(t) = \cos(2\pi f_c t).$$ The wireless channel will add a random phase, so the received signal is $$r(t) = \cos(2\pi f_c t + \phi)$$ (ignoring changes to the signal amplitude for the moment). Note that, since $\phi$ is random and changes over time, there is no real way for the receiver to determine it.

Now calculate the signal out of the quadrature branch of the receiver: $$ r_Q(t) = \text{LPF} \lbrace r(t) \sin(2\pi f_c t + \theta) \rbrace,$$ where $\theta$ is the phase of the receiver oscillator. If you do the math, you'll see that the energy of $r_Q(t)$ is not zero.


Low Pass Filter Artifacts I have created a Low Pass filter (I am trying to understand the topic of DSP, and it's hard.. But moving along with samples, books and a nice post on Stack Overflow).

I am saving the 'output' to a file, and I hear a constant clicking noise (very rapid, medium-pitch). The volume (intensity) of the 'clicking' gets louder as there is more intensity (in the passband of my filter).

I have been reading something about possible Artifacts when designing filters (FIR), but I am wondering if this is the constant clicking noise I hear ? It's not bad for my system, because I try to measure the intensity of the 0..100Hz band, and the clicking increases when the low sound increases. But I am just wondering if it's bad, or a possible symptom of a flaw in my filter ?

Thanks!

My Filter (of-course the getLPCoefficientsButterworth2Pole() is called first):

double ax[3];
double by[3];
double xv[3];
double yv[3];

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
double PI = M_PI;
double sqrt2 = sqrt(2);

double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
double QcWarp = tan(QcRaw); // Warp cutoff frequency

double gain = 1 / ( 1 + sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) );

by[2] = ( 1 - sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) ) * gain;
by[1] = ( 2 - 2 * 2 / ( QcWarp * QcWarp ) ) * gain;
by[0] = 1;

ax[0] = 1 * gain;
ax[1] = 2 * gain;
ax[2] = 1 * gain;
}

void process(MTAudioProcessingTapRef tap, CMItemCount numberFrames,
         MTAudioProcessingTapFlags flags, AudioBufferList *bufferListInOut,
         CMItemCount *numberFramesOut, MTAudioProcessingTapFlags *flagsOut)
{
OSStatus err = MTAudioProcessingTapGetSourceAudio(tap, numberFrames, bufferListInOut,
                                                  flagsOut, NULL, numberFramesOut);

if (err)
    NSLog(@"Error from GetSourceAudio: %ld", err);

float leftVolume, rightVolume;

NSLog(@"Test: %f", xv[1]);

for (CMItemCount i = 0; i < bufferListInOut->mNumberBuffers; i++) {
    AudioBuffer *pBuffer = &bufferListInOut->mBuffers[i];
    int cSamples = numberFrames * pBuffer->mNumberChannels;

    float *pData = (float *)pBuffer->mData;
    float rms = 0.0f;

    for (int j = 0; j < cSamples; j++) {

        xv[2] = xv[1];
        xv[1] = xv[0];

        xv[0] = pData[j];
        yv[2] = yv[1];
        yv[1] = yv[0];

        yv[0] =   (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                   - by[1] * yv[0]
                   - by[2] * yv[1]);

        pData[j] = 20 * yv[0];
        rms += yv[0] * yv[0];

    }

    if (cSamples > 0) {
        rms = sqrtf(rms / cSamples);
    }

    if (0 == i) {
        leftVolume = rms;
    }

    if (1 == i || (0 == i && 1 == bufferListInOut->mNumberBuffers)) {
        rightVolume = rms;
    }
}
}

From the audio clip posted on soundcloud, it looks like there's a discontinuity in the output every ~2060th sample. I suspect that you process your audio stream by blocks and that you do not save/restore the state variables of your filter between each call to the filtering function - so every time you process a new chunk of audio it is as if it was preceded by zeros - which is not the case obviously.

Related question.


ensure stability against coefficient quantization for a filter:

$$ H(z)=\frac{1+0.1z^{-1}}{1+0.1z^{-1}+0.9998z^{-2}} $$

Which precautions could be taken to ensure that the filter does not become unstable because of coefficient quantization? I've think that it would depend of the implementation chosen (I've read there are implementations pretty robust against quantization), but not really sure.  Thanks in advance.

your filter is in the form of a standard biquad with 5 independent coefficients ($b_2=0$ in your case):

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{ 1 + a_1 z^{-1} + a_2 z^{-2}} $$

your two poles are the value of $z$ that make the denominator equal to zero

$$ 1 + a_1 z_p^{-1} + a_2 z_p^{-2} = 0 $$

or (multiplying both sides by $z_p^2$),

$$ z_p^2 + a_1 z_p + a_2 = 0 \tag{1} $$

and stability depends entirely on whether your poles (after coefficient quantization) are contained in the unit circle:

$$ |z_p| < 1 $$

or

$$ |z_p|^2 < 1 $$

the solution to Eq. 1 is (quadratic formula):

$$ z_p = -\frac{a_1}{2} \pm \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} $$

and, if the poles are complex conjugate, then $\left(\frac{a_1}{2}\right)^2 < a_2$ and your complex conjugate poles are

$$ z_p = -\frac{a_1}{2} \pm j \sqrt{a_2 - \left(\frac{a_1}{2}\right)^2} $$

in the case of two real poles, when $\left(\frac{a_1}{2}\right)^2 \ge a_2$, you must guarantee (to ensure stability):

$$ -1 < -\frac{a_1}{2} \pm \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} < +1 $$

that means both of the following must be guaranteed (after quantization):

$$ -\frac{a_1}{2} + \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} < +1 $$

$$ -1 < -\frac{a_1}{2} - \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} $$

what is necessary (but not sufficient) is that $|a_1| < 2$ and what also must be true is (getting rid of the $\sqrt{\cdot}$)

$$ \left(\frac{a_1}{2}\right)^2 - a_2 < \left(1 + \frac{a_1}{2}\right)^2 = 1 + a_1 + \left(\frac{a_1}{2}\right)^2 $$

and

$$ \left(\frac{a_1}{2}\right)^2 - a_2 < \left(1 - \frac{a_1}{2}\right)^2  = 1 - a_1 + \left(\frac{a_1}{2}\right)^2 $$

or

$$ a_2 > -1 - a_1 $$

$$ a_2 > -1 + a_1 $$

so for real poles ($\left(\frac{a_1}{2}\right)^2 \ge a_2$) you must make sure, after quantization that

$$|a_1| < 2$$
and
$$-1 + |a_1| < a_2 \le \left(\frac{a_1}{2}\right)^2$$

for complex conjugate poles ($\left(\frac{a_1}{2}\right)^2 < a_2$), it's much, much easier.  i'll leave this for an exercise for the reader but you can easily show that (for complex conjugate):

$$ |z_p|^2 = a_2 $$

which means

$$ 0 < |z_p|^2 = a_2 < 1 $$

and you still will see that $|a_1| < 2$ is necessary because $|a_1| < 2 a_2 < 2$.  it turns out that $-\frac{a_1}{2}$ is the point on the $\mathfrak{Re}(z)$ axis that is midway between the two poles, whether they be real or complex conjugate.

now, if your digital filter is slewing from one stable complex pole case to another stable complex pole place, then $a_2$ is slewing from one value less than 1 to another value less than one and all filters in between will have $a_2$ less than one and will be stable, at least in a "static" sense.  rapid coefficient slewing can create instabilities in and of itself, but if the slewing stops (at a stable condition), any instability also stops.

so in both real and complex conjugate cases, to have stability in the above biquad you must have

$$ -2 < a_1 < 2 $$

and

$$ 0 < a_2 < 1 $$

for the complex conjugate case $a_1^2 < 4 a_2$

or

$$-1 + |a_1| < a_2 \le \left(\frac{a_1}{2}\right)^2 = \frac{|a_1|^2}{4}$$

for the real pair case $a_1^2 \ge 4 a_2$.

note that, in the real pole pair case, we know that $$-1 + |a_1| <  \frac{|a_1|^2}{4}$$ because $$-1 + |a_1| - \frac{|a_1|^2}{4} = -\frac{1}{4}(|a_1|-2)^2 < 0$$.  so there is always a gap for $a_2$ to fit into.


Is $x(t) = \cos t + \sin\left(\frac{1}{2}t\right)$ a periodic signal? Is $x(t) = \cos t + \sin\left(\frac{1}{2}t\right)$ a periodic signal?

The answer provided by the book is different from my answer. Book says its not a periodic signal. Can you guys tell me why is it not a periodic signal?

My answer:

  $\cos(t)$ is periodic as $2\pi f_1 = 1 \Rightarrow f_1=\frac{1}{2\pi} \Rightarrow T_1=2\pi$

  $\sin\left(\frac{1}{2}t\right)$ is also periodic as $2\pi f_2 = \frac{1}{2} \Rightarrow f_2=\frac{1}{4\pi} \Rightarrow T_2=4\pi$

  Therefore $\frac{T_1}{T_2} =  \frac{2\pi}{4\pi} = \frac{1}{2}$  is rational number

  Therefore the given $x(t)$ is a periodic signal.

As for every $t_0\in\mathbb{R}$ and $k\in\mathbb{Z}$
$$
\begin{eqnarray}
  &x(t_0+4k\pi)  &=\cos(t_0+4k\pi)+\sin(t_0/2+2k\pi)\\
  &              &=\cos(t_0)+\sin(t_0/2)\\
  &              &=x(t_0)
\end{eqnarray}
$$
you answer is correct: $x(t)$ is periodic.


Unexpected Result When Using Sinc Interpolation blue is how I tried to sinc interpolate. why would something like this happen?

Since Sinc based Interpolation requires you to know the data at any point it can not be done.

You might do a Truncated Sinc Interpolation.
The artifacts you're seeing can be caused by a kernel which is too short or the parameters aren't good.

In order to create a good Sinc kernel you need to know things about the Band Width of the signal and the Sampling Rate, did you took those into account?


Stability of the summation function Stability of the below function was asked.

$$
T(x[n]) = \sum_{k=n_0}^n x[n]
$$

I understand that $x[n]$ shouldn't exceed $M$ value. $x[n]\le M$. I mean, we are looking for if the input is bounded, what will the output be.

However, i didn't understand from where $|n-n_0|M$ came from in the below function.

$$
\left |T(x[n]) \right| \le \sum_{k=n_0}^n \left| x[n] \right| \le |n-n_0|M
$$

Well, we are summing over $k$ in the range $n_0$ to $n$. The number of elements in the sum is therefore $n-n_0$. We also know that the sum must be positive and that all terms in the sum ($x[n]$) have absolute value less than $M$.

To account for $n-n_0$ being negative, but our sum MUST be positive, we take its absolute value.


Is there a relationship between Average energy of transmission and power saving As mentioned  in many books, the average energy of transmitted symbol is given by:

$E_s = (2/3)(M - 1) $

where $M$ indicated the modulation order.
Suppose that we are using in the first scenario $SISO$ system, and we need to transmit 8bit/channel use, so the needed modulation is $256QAM$ which is equivalent to $E_s = 170$.  now suppose that we are using $MIMO$ system with two transmitters antennas and two receiver antennas, so to transmit the data rate of  8bit/channel use, we need to transmit 4bit/channel use in each antenna, we will need to use $16QAM$, so the needed $E_s = 10$ in each antennas, which is in total 20.

can we say, the power saving for transmitting 8bit/s/Hz in case of using MIMO system is 170/20  ?

Is that right? thank you

First of all, 8 b/s/Hz is not a bit rate, but rather a spectral efficiency.

Second, you don't transmit $\log_2(M)$ bits per second, but rather $\log_2(M)$ bits per channel use, which is roughly $T_s$ seconds, where $T_s$ is the symbol time. So, the bit rate is $\frac{\log_2(M)}{T_s}$.

If you adjust your wording to the above by saying 8 bits/channel use, then, yes, you are right, you save energy by using MIMO, but the saving is $170/20$ not the reciprocal.


How to get coordinates(x(time),y(frequency)) of spectral peaks from spectrogram? I have created this spectrogram from a wav file. Please have a look:

As it can be observed clearly that my somewhat spectral peaks are visible. I want to get the x = time, y= frequency of these spectral peaks in order to find temporal displacement between the peaks. How can I get these coordinates?

I am looking for a programmatic approach (using python,java etc).

Regards,

Khubaib

Under the hood, your spectrogram consists of multiple, consecutive FFT frames (of length N) and each FFT frame contains N/2 frequency bins (even if you have N frequency bins, the upper part is discarded).

An FFT frame denotes a temporal position from the beginning of a signal and this temporal position can be computed as

time_pos = frame_number *  N / sampling_rate; //frame_number runs from 0......N-1

As an example, if your signal is sampled at 8000 Hz and N is 1024, the 10th frame will be at 10 * 1024 / 8000  = 1.28 sec from the beginning of the signal.

Just like every FFT frame has a fixed length (N samples long), each frequency bin has a fixed width in Hz (frequency resolution). This frequency resolution is given by the sampling rate divided by the frame/fft length so

freq_res = sampling_rate / N;

In order to find the frequency for a particular bin, simply multiply the bin index by the frequency resolution.

freq_10bin = 10 * freq_res;

In order to get the coordinates of the spectral peaks, you should iterate through all the frames and find those bins that have a value above some threshold thus indicating a peak.
The frame indexes and bin indexes respectively are then used to compute the temporal and frequency coordinates for every peak (and the corresponding displacements).


How to select the length of a smoothing moving average filter? I have a time series signal I need to smoothen to get rid of some noise.

I decided to apply a moving average filter to this signal. I know the choice of the length of the MA is critical.

How should I choose the length?

The magnitude frequency response $|Y(\omega)|$ of moving average of length $N$ can be approximated by:

$$|Y(\omega)| \approx |\hat Y(\omega)| = \left|\frac{2\sin(\frac{N\omega}{2})}{N\omega}\right|$$

Figure 1. Horizontal axis: $N\omega$, vertical axis: $|\hat Y(\omega)|$.

The frequency $\omega$ is in radians, meaning that $\omega=2\pi$ corresponds to the sampling frequency. You should choose $N$ such that the noise is attenuated as much as possible while keeping a flat enough frequency response over the signal's bandwidth so that the signal's waveform is not distorted too much. What is flat enough depends on what you are going to do with the signal. A good first bet is to have the -3 dB cutoff frequency at the upper bandlimit of your signal so that no frequency of the signal is attenuated by more than 3 dB. The -3 dB cutoff frequency $\omega_c$ of moving average can be approximated by:

$$\omega_c \approx \hat\omega_c = \frac{2.783}{N}$$

These approximations were derived here.


How to calibrate three orthogonal cameras？ I am trying to find out intrinsic parameters(focal length ,principle point,and lens distortion coefficients)as well as the extrinsic parameters(rotation and translation of each camera relative to the reference frame) of three cameras in the flowing picture.

They are placed at nearly orthogonal angles around a 12x12x12 inch working space. There are four steps:

capture images of chess board by two cameras at the same time
use toolbox calib_gui(matlab) to calibrate two cameras separately and get Calib_Results_left.mat and Calib_Results_right.mat.
put these two results into stereo_calib(matlab) and run Calibration
repeat step 1~3 on another pair of cameras

However, I was stuck at step 3! None of pair images was found inconsistent. The errors were shown in the following form:

Disabling view 2 - Reason: the left and right images are found inconsistent (try help calib_stereo for more information)

I checked the order of images but didn't find problems. I doubt that my application is not binocular calibration, because intersection angle between two cameras is too big
for a stereo-camera!

I just want to find out rotation and translation of each camera relative to the reference frame. Does anybody have any ideas?

Use an asymmetric calibration plate. Such as the new circle grid in OpenCV. Check here. You do not need to find the correspondences in the gui, but rather run the calibration directly. Also, do not forget to initialize the optimization problem using the median of the obtained poses.

To calibrate multiple view setups, I would also recommend Multi Camera Self Calibration tool. It uses a simple laser dot and is much easier to calibrate, especially when the number of cameras increase. Of course, pay attention to the synchronization.


Machine learning with only positive labels Suppose I have a binary classification problem with 10 features and about 1000 samples. In the training set, most of my data is unlabeled (75%). The rest of the data is labeled but contains only positive labels.

In the test set, I have both negative and positive labels. How should I approach this classification problem?

I would use a novelty detection approach: Use SVMs (one-class) to find a hyperplane around the existing positive samples. Alternatively, you could use GMMs to fit multiple hyper-ellipsoids to enclose the positive examples. Then given a test image, for the case of SVMs, you check whether this falls within the hyperplane or not. For GMMs, you check if it is enclosed in the hyper-ellipsoids. They are both proven to work well in practice.

If you also have some unlabled data in your training set, I would certainly adapt a variant of transfer learning. Maybe you would be able to automatically label the unlabeled data, based on the already learnt samples.