How to perform a phase recovery for higher order modulation schemes My objective is to demodulate the signals with different modulation orders and schemes. In particular, I want to recover the phase first.
Let us say, my signal looks like in the constellation below
.
The idea is to position the constellation in such a way that the largest excursions are at 45 degrees angle. In other words, IQ diagram must be shifted by some angle. This angle is determined by the average shift of all points from the reference axis (real or imaginary). The idea itself is coming from this post, where BPSK phase recovery is described.
My problem is that I still do not understand, how do we perform the phase recovery for the higher order modulation schemes. The idea of simply mirroring the constellation to one quadrant (as it works in the case of BPSK) does not seem correct to me and also does not work well. So how we can do it?
Thanks!

This post provides many details on implementing carrier recovery for a higher order QAM system (or even QPSK).
High modulation index PSK - carrier recovery
Using a PI loop filter where in addition to a simple accumulator, a direct (proportional) path is also summed in results in a 2nd order Type 2 tracking loop as required to track phase offsets to zero error under static frequency offset conditions. Symbol-aided synchronization (where a known pilot is used) is recommended for higher order QAM to minimize acquisition time. The coefficients of the loop filter are set to stabilize the carrier recovery loop and also set the appropriate loop BW.
See this post for further considerations on setting the loop BW properly:
Loop bandwidth for symbol timing recovery


Stability of a logarithmic function and system A system such that

$$y[n] = \log_{10}(|x[n]|) $$

was given as a stable one. I am not able to understand it though. When $x[n] = 0$ (which is bounded), the output $y[n]$ tends to get unbounded. So how is this system stable?

As correctly pointed out in a comment, this system is definitely not BIBO-stable. The output signal becomes arbitrarily large (in magnitude) when $x(t)$ approaches zero. So a bounded input signal can result in an unbounded output signal, and, consequently, the system is unstable.


Synchronicity of two signals I am trying to determine the synchronicity between two oscillating signals. I was attempting to do this by making use of the cross correlation between the two signals and looking at the maximum value. Is this a correct approach?

If so, which would be the better approach, to normalize the signals between 0 and 1 or to normalize the cross correlation?

I don't really mind about the difference in amplitude but I would like to see levels of synchronicity. Originally I thought to do cross correlation and look at the max of the graph however, I am not sure that this is the correct approach. I have also been looking at coherence of the spectra

Using cross-correlation and finding the peak is one way of finding the time delay between two signals.
This assumes that the two signals are periodic and have the same fundamental frequency.

When the two signals have close, but not equal fundamental frequency, the lag position of the maximum cross correlation value will change in time.  If the lag of the maximum does not change (and some thresholds are satisfied so you're not just cross-correlating two very weak signals with the same 60 Hz hum in the background), the two signals are synchronous.  If it changes slowly, that indicates the degree that the two signals have different frequencies.  If the lag requires T seconds for the lag to slide one period (which would be the distance between lags of maximum correlation) in one second, the two signals differ in frequency by 1 Hz.


How to implement Discrete time FIR Hilbert transformer filter? I want to implement FIR Hilbert transformer filter in Matlab.

From your comment it became clear that you're actually interested in the design of an FIR Hilbert transformer. There are many ways to design a Hilbert transformer, but for you one of the most straightforward ways is probably to use the Matlab function firpm.m with 'ftype' set to 'h' (for Hilbert transformer). This will design an FIR Hilbert transformer with an equiripple behavior (i.e., with the maximum approximation error being minimized). Alternatively, you could use firls.m for a least squares design. If you use firls.m you need to set the filter type 'ftype' to 'hilbert'.


Confused about the definition of point to point ,this point means single antenna or single device? What is the definition of point to point ? does the "point" mean single device or single antenna?  I mean,can one transmitter with multiple antennas and one receiver with multiple antennas be called point to point ?or we should only call it MIMO?

If there are 2 transmitters,t1 & t2,and 2 receivers,r1 & r2,and t1 will send the message to r1,but r2 will also receive the RF of t1 ;t2 will send the message to r2,but r1 will also receive the RF of t2.So in this situation,is this system also point to point? or point to multi-point

If a receiver with multiple antennas ,and there are K energy harvesters, 4 antenna for each energy harvester,and 1 receiver with single antennas.What system can we call it?MISO?or MIMO?

Point to point typically refers to single devices transmitting to each other, regardless of how many antennas they have. The important thing is that in point to point systems, there is no interference from other transmissions, you only have signal and noise.

If transmitter and receiver have multiple antennas, you would typically call it MIMO, to specify even further, the terms single user MIMO (SU-MIMO) or point to point MIMO (P2P-MIMO) are not uncommon.

If only one antenna is used for data and the others for harvesting energy, it depends on your viewpoint. If it makes sense to treat these antennas separately, i.e., ignoring the energy harvesting antennas when designing the communication link, there might be a point to call it MISO. I would assume that one also optimizes the energy transfer so that all antennas are looked at jointly. This would advocate for calling it MIMO.

*edit: Regarding the scenario you edited in later: As Blackmath correctly said, this is a classical interference channel. So it's not a point to point setup. It's also not Point to Multipoint, in which the same information is broadcast to several users. Not so in your setup as you have different intended pieces of information for your two receivers. This is an interference channel.


Linear Chirp parameters from a segment of the signal Edit:" please see the attached figure. The blue part is all I have. Also, you can neglect the noise term. Assume the signal is deterministic. Typical values of the unknown parameters $\alpha$ and $\omega_0$ are in the the range [0.25,5]. Note that all the parameters are unknown $A \omega_0, \alpha, \phi_0$. The setup hints at curve fitting which I did but don't prefer."

I have a small part from a chirp signal

$$s(t)= A \cos\left(\omega_0 t+ \frac\alpha2 t^2 + \phi_0\right)+n(t)$$

with very low start frequency $\omega_0$  and chirp rate $\alpha$.

The signal is sampled at high sampling frequency $f_s$. But the available time domain signal is short, i.e. barely a complete cycle. The amplitude of the signal and the phase shift are also unknown. Are there ways to figure out the chirp parameters other than curve fitting or spectrogram?

Amplitude $A$ is trivial. Square and low-pass filter/average; then take the square root.

Without noise (or with good SNR), you could simply calculate a derivative:

\begin{align}\frac{\partial\, s\}{\partial t} = -A\sin\left(\omega_0t + \frac\alpha2t^2 +\phi_0\right)\cdot\left(\omega_o+ \alpha t\right)
\end{align}

and low-pass filter the result

$$\text{LPF}\left\{\frac{\partial\, s\}{\partial t}\right\} = \omega_o + \alpha t$$

From the resulting linear function, $\alpha$ can directly be extracted as the slope, and knowing the slope, $\omega_0$ is just the constant summand.

With $A$, $\omega_0$ and $\alpha$ known, you can simply compare at any fixed $t$ to find $\phi_0$.


Understanding BFSK Modulation (MATLAB/Simulink) Could someone explain to me BFSK Modulation on the Simulink model below?

The simulation above simulates a binary FSK system with a correlation (or matched filter) receiver.

I have to understand the general idea about it.

Could someone explain this?
Moreover I cannot understand the meaning of correlator and tone frequency fn the above simulator?

So, ignoring noise, the input into the receiver multiplier coming from the channel is either $\cos(2\pi f_0t)$ or $\cos(2\pi f_1t)$. The other multiplier input is $\cos(2\pi f_0t)-\cos(2\pi f_1t)$. All signals are defined from 0 to $T$. Let's assume the transmitted signal has frequency $f_0$. Then you have an integrator:

\begin{align}
&\int_0^T \cos(2\pi f_0t) \left[\cos(2\pi f_0t)-\cos(2\pi f_1t)\right]dt \\
&=\int_0^T \cos(2\pi f_0t)\cos(2\pi f_0t) dt - \int_0^T \cos(2\pi f_0t)\cos(2\pi f_1t)dt.
\end{align}

The frequencies $f_0$ and $f_1$ are chosen so that the cosines are orthogonal for $0<t<T$; that means that
$$\int_0^T \cos(2\pi f_0t)\cos(2\pi f_1t)dt =0.$$
The integrate and dump block outputs the result of the definite integral at $t=T$. You can see that the result is positive if frequency $f_0$ was transmitted and negative if $f_1$ was transmitted.

The sign block is a quantizer, and the lookup table converts the numbers back to bits 0 and 1.


Can error in digital signal be zero after filtering? Can error in digital signal be zero after filtering? I have looked up online to know about this but I could not find a solution.

The discrete signal is a sum of two sines. The noise is an additive uniform noise, always positive (rand lies in open interval $(0,1)$). The latter  is the crucial point with the median. Suppose you have a frame of three signals samples: $x_0$, $x_1$, $x_2$, and three noise realizations:  $n_0$, $n_1$, $n_2$ in $[0,1]$. So you are checking if:

$$\textrm{med}\{x_0+n_0,x_1+n_1,x_2+n_2\}-x_1$$

can be zero for each window frame. You have noticed from the plots that the noisy signal is always above the clean one. So the  median of $\{x_0+n_0,x_1+n_1,x_2+n_2\}$  is always higher than the  median of $\{x_0,x_1,x_2\}$.

But could it be the same as the initial data (putting signal borders away)? Yes. Take  $\{x_0,x_1,x_2\}=\{0,0.5,0.2\}$ and add the random realization $\{0.5,0.1,0.3\}$. Then the median of the noisy signal is $0.5$, the exact value of the center sample $x_1$.

So for some well-chosen signals and noise, their are counter-examples, a median filter can in cases provides a zero error. However, and I admit I cannot give a proof at this time, it is probably very unlikely that there is a zero error for  long enough generic signals.

Note to DSP folks: I feel like I am missing a very evident proof, or such a question is not adapted at this level of signal processing.


Spectral averaging without phase information In the context of signal averaging to reduce uncorrelated noise:

Is it possible to perform signal averaging from just the magnitude spectrum of certain signals? Or that would only work if we know also the phase information of those signals. In that case, what would be the requirement for the signals to average, regarding phase information? I guess they both have to be somehow "related" to the same "starting point".

The short and most common answer is that you can average without phase but there are cases where you can average phase as well.

Cross correlation between
2 sigmals where one has a delayed component of the other that is of a "long" duration can exploit complex averaging.

Spectral analysis using the periodogram averages magnitude or magnitude squared.

There are numerous examples of where you do or don't use complex averaging.  It depends on the problem being solved and is typically going to be relatively obvious if you should or shouldn't.  If you expect the phase to evolve in a predictable way over your observations, you can use that to your advantage.


Discrete Wavelet Transform output: coefficients or FIR-filtered signals? I have a theoretical question about the calculation of the Discrete Wavelet Transform, using MATLAB specifically.

According to this video-tutorial on the MATLAB lagorithm: https://it.mathworks.com/videos/understanding-wavelets-part-2-types-of-wavelet-transforms-121281.html

the DWT algorithm applies a low-pass and a high-pass filter to a signal to obtain a low frequency signal and a high frequency signal, using appropriate filters. And this makes sense to me, because the output of the application of a filter, such as a FIR filter, (i.e using the FILTER function in MATLAB) is another signal.

However, looking at the documentation of the function DWT or WAVEDEC in Maltab, it looks like the output of the j-th step of the DWT is NOT another signal, rather the cAj and cDj coefficients, from which the low-frequency and high-frequency signals can be then reconstructed.

So my question is: what is the filter that is applied to the signal? Why does it provide coefficients and not another signal? What is the difference between the application of the DWT and that of a FIR filter?

In other words, it seems that I apply a filter (with some coefficients) to find other coefficients, rather than another signal ... this is not 100% clear to me.
Maybe I misunderstood some basic concepts, can you please clarify?

In one level DWT, each output of the low-pass or a high-pass can indeed be considered as signals. Thus each of those signals are subsampled by a factor of 2, and the same two-filter-subsampling is iterated on the low-pass output, several times (wavelet decomposition) at $L$ levels. Each final output of the different branches could still individually be considered as "signals", but they only make sense together, with respect to the input signal. Coefficients is a traditional name for them. Honestly, I don't know where it appeared originally.

Each represents information cast at a certain scale, and they don't really belong to the same domain. They are signals projected onto a nested wavelet subspace from a multiresolution analysis. They altogether belong to a global transformation. In other words: the values observed are signals with a heavy wavelet footprint (especially because of the subsampling). Analyzing them alone (without taking into account the wavelet projection) is theoretically dangerous, and even called "wavelet sin", or "wavelet crime", see for instance: Wavelet transform: How to compute the initial coefficients when only samples are available? (Strang and Nguyen, Wavelets and Filter banks, 1996, pages 232 sq.). However, people often may sin with practical efficiency, and for orthogonal non-redundant transformation, it often does not go too bad.

Anyway, due to the invertibility of the DWT scheme, you would have "more genuine" signals by: keeping some coefficients (setting the others to zero), especially from each "subband coefficient set alone", and performs the (unique) inverse DWT for each of the  $L+1$ sets of subband coefficients. From those $L+1$ inverted, you get $L+1$ reconstructed signals in the original domains, whose sum yields the original signal (linearity of the wavelet).

The "wavelet crime" often becomes more stringent when one uses union of wavelet bases, wavelet frames, etc. As for the MODWT, I would not claim (yet) for sure that I fully understand the concept. Its setting sis akin to undecimated, shift-invariant cycle-spinning or redundant discrete wavelets. While the DWT generates $N$ coefficients for a $N$-sample signal, whatever the levels, the MODWT generates about $(L+1)N$ coefficients. As far as I get it, MODWTMRA does a job similar to the one described above for the DWT: project the redundant wavelet coefficients back to th original domain, so that their sum yields the original signal back. Side note: for redundant transforms, the inverse is not unique.


The meaning of $k \delta (t)$ Reading the signal literature I often come across the expression $k \delta (t)$ where $k$ is a constant. I presume this is a notation to suggest we are referring to the area or strength of the Dirac function since the multiplication of a constant, $k$, by $\delta (t)$ seems to make no sense.

My question is what does $k \delta(t)$ mean? Am I correct in assuming that on its own it doesn't represent $k$ multiplied by $\delta(t)$?

You are correct. In your example, $k$ would refer to the "area underneath" the impulse. Mathematically speaking, the Dirac delta isn't a function in the typical sense of the term. Instead, it is more of a distribution, characterized by the fact that when integrated across any interval that contains $t=0$, the result is unity. That is:

$$
\int_{-\epsilon}^{\epsilon} \delta(t) dt = 1 \ \forall\ \epsilon > 0
$$

The distribution is typically defined as follows:

$$
\delta(t) = \begin{cases} \infty,\ t = 0 \\ 0, \text{ otherwise}\end{cases}
$$

Multiplication by a constant $k$ obviously wouldn't change this definition at all. So instead, the scaling merely changes the area underneath the distribution:

$$
\int_{-\epsilon}^{\epsilon} k\delta(t) dt = k \ \forall\ \epsilon > 0
$$

Extending this concept of multiplication by a constant $k$ to multiplication by a function $f(t)$ yields the following:

$$
\int_{-\epsilon}^{\epsilon} f(t) \delta(t) dt = f(0) \ \forall\ \epsilon > 0
$$

or more generally:

$$
\int_{T-\epsilon}^{T+\epsilon} f(t) \delta(t-T) dt = f(T) \ \forall\ \epsilon > 0
$$

which is known as the sifting property of the Dirac delta and is used extensively in linear systems theory.


What's the advantage of adaptive IIR filter against FIR? Adaptive IIR filters is not straightforward, and may be unstable. Many people say that adaptive IIR filters use less coefficients than FIR filters. What I'm curious about is how many coefficients can IIR save?

I tried to use adaptive IIR filters to estimate transfer function of a 32-order FIR filter. Assume the IIR filter has $M+N+1$ coefficients: $a_1, a_2, ..., a_M, b_0, b_1, ...b_N$. I found the estimation result is acceptable only when $M+N+1 \ge 30$, i.e. only 2 coefficients can be saved.

In actual projects, for example, a 50 MHz FPGA, a 32-order FIR will produce about $(32 / 50 ~{M}) / 2 = 0.32 ~{\mu s}$ delay, so

What will happen for IIR?
Can adaptive IIR filters really reduce number of coefficients and reduce signal processing time delay?

These are the key differences between FIR and IIR filters, regarding the feature you wish to control are the following:

$$
\begin{array}{c|lcr}
\text{Feature} & \text{IIR} & \text{FIR} \\
\hline
\text{Implementation} & \text{Poles & Zeros} & \text{Zeros Only} \\
\text{States} & \text{Yes} & \text{No} \\
\text{Phase Delay} & \text{*} & \text{Half Integer} \\
\text{Stability} & \text{*} & \text{Always} \\
\text{Ripple} & \text{Yes} & \text{*} \\
\text{Cut-Off} & \text{Yes} & \text{*} \\
\end{array}
$$

The * indicates the feature can be controlled, by adding orders in most cases.

The standard definitions of FIR and IIR filters are:

FIR:

$$H(z)=b_0z^0+...+b_nz^n$$
$$y(t)=b_0u(t)+...+b_nu(t-n)$$

IIR:

$$H(z)=\frac{b_0+b_1z^1+...+b_nz^n}{1+a_1z^1+...+a_nz^n}$$
$$y(t)=b_0u(t)+...+b_nu(t-n)-a_1y(t-1)-...-a_ny(t-n)$$

$u$ is the input, $y$ is the output, $x$ is the states (below), $t$ is the time, scaled by a sampling time $dt$, $n$ is the number of orders of the filter. Each filter has $n$ size coefficient vectors, plus constant direct-output term $b_0$ (optional), and $a_0$=1. For simplicity assume $\sum{b_i}=1$ and $\sum{a_i}=1$, though this is not required anywhere.

Implementation. By definition, FIR include zeros only, leading to a linear system in the history vector for $u$: $[u(t-1)...u(t-n)]$.

IIR include both poles and zeros, also leading to a linear system in the history vector not only for $u$, but for $y$ too. Because of this, by one side IIR can be unstable; but by other side, they can be designed to have smooth ripple and sharp cutoffs with a minor number of orders.

States. FIR are static systems in the history vectors, meaning the filter is not dynamical, do not have states, is not recursive, no feedback.
IIR are dynamical systems in the history vectors, meaning the filters have states, is recursive, has feedback, hence have "memory" from past inputs & outputs.

Phase Delay. The Phase Delay $\tau_\phi$

$$y(t)=y_0(t-\tau_t) sin(\omega(t-\tau_\phi)+\theta)$$

can be easily controlled in FIR implementations.
If $b_k=b_{n-k}$,$k=0...n$, the phase delay is constant, equal to $n/2$ (the center of the FIR coefficients shape, its impulse response), equal to the group delay, and thus the filter become linear phase, with phase equal to $\omega\tau_phi$.

Because IIR have infinite impulse response, they can be minimum phase instead of linear phase, though the phase achieved can be much less than the phase of a FIR for the same number of orders.

Stability. FIR are always stable, IIR can be designed to be stable, if stability is required.

Ripple. IIR can be designed to be flat-ripple both in pass-band|stop-band|both (butterworth|chebyshev|elliptic), FIR requires a major (tending to "infinite") number of orders for equate this property.

Cut-Off. IIR can be designed to have a sharp cut-off or narrow transition bands, FIR requires a major (tending to "infinite") number of orders for equate this property.

Related Articles:

https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-341-discrete-time-signal-processing-fall-2005/lecture-notes/lec08.pdf
https://www.quora.com/Why-are-FIR-filters-preferred-over-IIR-filters
http://iowahills.com/A8FirIirDifferences.html
http://forums.prosoundweb.com/index.php?topic=2045.0
http://www.vyssotski.ch/BasicsOfInstrumentation/SpikeSorting/Design_of_FIR_Filters.pdf


How to filter objects on basis of orientation of regionprops I am trying to detect straight vertical lines in an image. I used vertical edge detection and then I tried to check which objects has the orientation > 80 & < 100.

Input Image:

Code:

% Vertical Edge Detection

v = edge(img,'sobel','vertical');
v = bwlabel(v);
stats = regionprops(v, 'All');

% ids = find([stats.Orientation] > 80 & [stats.Orientation] < 100);
% v = ismember(v,ids);

figure,imshow(v);title('Vertical');

Problem:

The problem is that most of the objects are having some weird orientation (some -ve value e.g. -53.xyz)

Output:

when I try to filter the objects, it gives me a full black image as an output.

Someone please explain this behavior. How can I get only straight vertical lines?

I don't have MATLAB, but here are a few ideas:

Since you use regionprops, you could use the Orientation parameter - which is pretty much exactly what you want to have.

Another idea

Create an matrix of the same size as your image. The matrix has entries that give a linear variation in the perpendicular orientation that you desire (it basically is a gradient field in the perpendicular orientation, e.g. with values ranging from 0 to 1).
Use bwlabel to label all your candidate lines
Loop over all labels and create a binary mask for each label. Multiply the mask with the gradient matrix.
Calculate the variance inside the mask. Since the gradient is perpendicular to your direction, the variance should be 0 for a perfect line. For a not-so-perfect-line, there is a certain variation. You could use

an empirical threshold
check the math for an more elaborate threshold, because one can certainly find a function $\sigma^2(\alpha)$, with $\alpha$ being the angle of the region of interest in the mask.


Artifact in Bicubic interpolation I wrote an algorithm to do bicubic interpolation of an image.
I used the method desribed in the wikipedia page.
On simple images, the result looks good, but on more complex ones, I got strange artifacts in non-smooth zones.

What could be a probable source of those artifacts, and what can I try to get rid of them ?

My code is on github, but I checked it several times, so I don't think the problem comes from it but rather from the simple implementation suggested in wikipedia

I'm using this particular equation (3rd-order Hermite polynomial interpolation) to do the interpolation (in one dimension):

$$x(n+t)=\frac{1}{2}
\begin{bmatrix}
1 & t & t^2 & t^3 \\
\end{bmatrix}
\begin{bmatrix}
0 & 2 & 0 & 0 \\
-1 & 0 & 1 & 0 \\
2 & -5 & 4 & -1 \\
-1 & 3 & -3 & 1 \\
\end{bmatrix}
\begin{bmatrix}
x[n-1] \\
x[n] \\
x[n+1] \\
x[n+2] \\
\end{bmatrix}$$

where $x[n-1]$, ... $x[n+2]$, .. are the values of the surrounding pixels (in one dimension), and $t$ is the fractional coordinate (where $0 \le t < 1$) of the interpolated pixel between the adjacent pixels $x[n]$, $x[n+1]$. $n$ is an integer and $x(n)=x[n]$ for any integer $n$.

The expected output:

My result, look closely at the bird head,a lot of white pixels have appeared.

---EDIT---
To show that my code does produce coherent results, here is the example of my scaling up on simple image with 16 red pixels.

And here is the ouput of the library (which also use 3rd-order Hermite polynomial interpolation):

And here is my results:

In this case the difference are IMO totally acceptable, but it shows that my interpolation does something legit.

What you're seeing is acutance, as described in the wikipedia article on bicubic interpolation.

However, for downsampling an image, cubic interpolation is perhaps not the ideal choice, but makes more sense for upsampling.

Note that if you're downsampling by more than a facor of 4, there will be many pixels in the original image that are simply never used in the downsampled version.

The easiest way to get around this is to downsample by iteratively taking the mean value of 4 neighbouring pixels, until you're within a factor of 2 of the desired image size. Then you can do the last part by cubic interpolation.

If you're want high quality resampling, try taking a look at the ImageMagik documentation for inspiration. In particular they recommend using windowed filters for downsampling, such as Lanczos or Sinc.


Pole Zero plot given a Transfer function I've been looking at how to plot zeros/poles based on a transfer function. I found a couple of Tutorials online.

In the first youtube tutorial, the author brilliantly explains how to plot the zeros/poles.

In the second tutorial, the author explains how to obtain the magnitude characteristics of the frequency response, which is also fantastic. In the second tutorial, under the Example in "3.2.3 Transfer function of discrete-time systems", we have the following transfer function

Based on the transfer function, the poles and zeros can be defined as,

a = [1 -2.2343 1.8758 -0.5713]
b = [0.0088 0.0263 0.0263 0.0088]

This is where my confusion starts. based on the first tutorial, i'll have to plot all the zeros/poles along the x-axis (Or am I mistaken?).

But based on the MATLAB command to plot pole and zeros,

zplane(a,b)

I get this plot

The plot for the poles and zeros are scattered all over. How can i plot the poles and zeros manually in the Z-plane given my poles and zeros and obtain a similar output to MATLAB ?

Thanks for your help.

What you have are not the poles and zeros, but simply the filter coefficients, i.e., the coefficients of the numerator and denominator polynomials.

The poles are the roots of the denominator polynomial, and the zeros are the roots of the numerator polynomial. In Matlab they can be found by using the roots command:

p = roots(a);
z = roots(b);

Note that in general, poles and zeros are complex numbers, that's why they are plotted in the complex plane.

Just a remark: you used the zplane command with numerator and denominator interchanged, that's why the plot shows the zeros as crosses on the unit circle, and the poles as 'O's inside the circle. The correct call is zplane(b,a).


I wish to create a GNU Radio program in C++ that doesn't need GRC or Python, is there a reference or template for this? I want to create a functionally equivalent C++ file to the Python top_block.py file that is generated by GRC when executing a (blank example) flow chart.

I know that C++ can be used to define blocks specifically, and so I would assume that by creating a C++ file and including and using the methods from defined blocks through the relevant header files would allow for, essentially, the creation of GNURadio programs without the use of GRC or Python.

Has anyone created similar programs to this? Is there any chance I can make GRC generate such a file, or can someone point me to reference programs / documentation?

At it's core my question is how would I make the following Python file in C++?

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Top Block
# Generated: Wed Oct 16 17:45:58 2019
##################################################

if __name__ == '__main__':
    import ctypes
    import sys
    if sys.platform.startswith('linux'):
        try:
            x11 = ctypes.cdll.LoadLibrary('libX11.so')
            x11.XInitThreads()
        except:
            print "Warning: failed to XInitThreads()"

from PyQt4 import Qt
from gnuradio import eng_notation
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser
import sys

class top_block(gr.top_block, Qt.QWidget):

    def __init__(self):
        gr.top_block.__init__(self, "Top Block")
        Qt.QWidget.__init__(self)
        self.setWindowTitle("Top Block")
        try:
            self.setWindowIcon(Qt.QIcon.fromTheme('gnuradio-grc'))
        except:
            pass
        self.top_scroll_layout = Qt.QVBoxLayout()
        self.setLayout(self.top_scroll_layout)
        self.top_scroll = Qt.QScrollArea()
        self.top_scroll.setFrameStyle(Qt.QFrame.NoFrame)
        self.top_scroll_layout.addWidget(self.top_scroll)
        self.top_scroll.setWidgetResizable(True)
        self.top_widget = Qt.QWidget()
        self.top_scroll.setWidget(self.top_widget)
        self.top_layout = Qt.QVBoxLayout(self.top_widget)
        self.top_grid_layout = Qt.QGridLayout()
        self.top_layout.addLayout(self.top_grid_layout)

        self.settings = Qt.QSettings("GNU Radio", "top_block")
        self.restoreGeometry(self.settings.value("geometry").toByteArray())

        ##################################################
        # Variables
        ##################################################
        self.samp_rate = samp_rate = 32000

        ##################################################
        # Blocks
        ##################################################

    def closeEvent(self, event):
        self.settings = Qt.QSettings("GNU Radio", "top_block")
        self.settings.setValue("geometry", self.saveGeometry())
        event.accept()

    def get_samp_rate(self):
        return self.samp_rate

    def set_samp_rate(self, samp_rate):
        self.samp_rate = samp_rate

def main(top_block_cls=top_block, options=None):

    from distutils.version import StrictVersion
    if StrictVersion(Qt.qVersion()) >= StrictVersion("4.5.0"):
        style = gr.prefs().get_string('qtgui', 'style', 'raster')
        Qt.QApplication.setGraphicsSystem(style)
    qapp = Qt.QApplication(sys.argv)

    tb = top_block_cls()
    tb.start()
    tb.show()

    def quitting():
        tb.stop()
        tb.wait()
    qapp.connect(qapp, Qt.SIGNAL("aboutToQuit()"), quitting)
    qapp.exec_()

if __name__ == '__main__':
    main()

Apologies as I am unsure if this question is better suited here or on stack overflow, feel free to mark it for migration if it is not appropriate here!

you can create flowgraph in c++ like gqrx and i think GNU Radio Manual and C++ API Reference documents will help you. This is gnuradio c++ document for top_block if you could not use it let me know to write a sample for you. These documents are for gnuradio version 3.7 and in gnuradio 3.8 there is c++ code generation and gnuradio blocks are porting to c++ and you can track them in this link.
i hope it helps.


How is WSNR (weighted signal to noise ratio) defined? I'm studying about comparison between images to determine quality. I've learned about MSE, SNR and PSNR, and now I'm trying to understand WSNR, which I assume is similar to SNR but with weights.
I have a Python code to calculate the WSNR between two images, but I don't fully understand it and I don't know what is the mathematical definition of this measure. I can't find any paper or book that explain how it is defined and calculated.
Question:
How is WSNR (weighted signal to noise ratio) defined?
Python Code:
import numpy as np
from scipy.ndimage.filters import gaussian_filter as __gaussian_filter
from scipy.ndimage.filters import convolve as __convolve
from scipy.ndimage.filters import correlate as __correlate
from scipy.fftpack import fftshift as __fftshift

def wsnr(reference, query):
"""Computes the Weighted Signal to Noise Ratio (WSNR) metric.

    value = wsnr(reference, query)

    inputs
    ----------
    reference: original image data.
    query    : modified image data to be compared.

    output
    ----------
    value    : wsnr value
    """
    def __genetate_meshgrid(x, y):
        f = lambda u: u / 2 + 0.5 - 1
        (H, W) = map(f, (x, y))
        return (H, W)

    def __create_complex_planes(x, y):
        (H, W) = __genetate_meshgrid(x, y)
        (xplane, yplane) = np.mgrid[-H:H + 1, -W:W + 1]
        return (xplane, yplane)

    def __get_evaluated_contrast_sensivity(plane):
        w = 0.7
        angle = np.angle(plane)
        return ((1.0 - w) / 2.0) * np.cos(4.0 * angle) + (1.0 + w) / 2.0

    def __get_radial_frequency(x, y):
        (xplane, yplane) = __create_complex_planes(x, y)
        nfreq = 60
        plane = (xplane + 1.0j * yplane) / x * 2.0 * nfreq
        s = __get_evaluated_contrast_sensivity(plane)
        radfreq = abs(plane) / s
        return radfreq

    def __generate_CSF(radfreq):
        a = -((0.114 * radfreq) ** 1.1)
        csf = 2.6 * (0.0192 + 0.114 * radfreq) * np.exp(a)
        f = radfreq < 7.8909
        csf[f] = 0.9809
        return csf

    def __weighted_fft_domain(ref, quer, csf):
        err = ref.astype('double') - quer.astype('double')
        err_wt = __fftshift(np.fft.fft2(err)) * csf
        im = np.fft.fft2(ref)
        return (err, err_wt, im)

    def __get_weighted_error_power(err_wt):
        return (err_wt * np.conj(err_wt)).sum()

    def __get_signal_power(im):
        return (im * np.conj(im)).sum()

    def __get_ratio(mss, mse):
        if mse != 0:
            ratio = 10.0 * np.log10(mss / mse)
        else:
            ratio = float("inf")
        return np.real(ratio)

    if not len(reference.shape) < 3:
        reference = __convert_to_luminance(reference)
        query = __convert_to_luminance(query)
    size = reference.shape
    (x, y) = (size[0], size[1])
    radfreq = __get_radial_frequency(x, y)
    csf = __generate_CSF(radfreq)
    (err, err_wt, im) = __weighted_fft_domain(reference, query, csf)
    mse = __get_weighted_error_power(err_wt)
    mss = __get_signal_power(im)
    ratio = __get_ratio(mss, mse)
    return ratio

There might be several versions, yet the most probable is as follows. An SNR is an energy-dependent measure. In the time or spatial domain, it is agnostic to location: all samples have the same "weight" (uniform weighting). As an energetic measure, it can be computed in the Fourier domain as well (or any orthogonal transformation).

Then, from the sensory system, we know that in hearing, or vision, not all frequencies have the same perceptual weight: a human hearing a loud 35 kHz tone won't probably not notice it, while a bat or a dog will notice. Hence, the relative importance of a sound should be weighted by the perceptual sensitivity of the subject (bat, dog, human). This yields a different weighting across different frequencies. One shall not care about a 35 kHz tone for a human ear, and give it a little weight, compared to the [1-2] kHz band. Same happens for vision. Such weightings are used for instance in JPEG compression.

So instead of computing, in the frequency domain some:

$$-10 \log \sum_k \frac{\left(\hat{s_k}-\hat{n_k}\right)^2}{\hat{s_k}^2}$$
with $k$ an index of frequency bins, you use a predefined weighting for each $k$:
$$-10 \log \sum_k w_k\frac{\left(\hat{s_k}-\hat{n_k}\right)^2}{\hat{s_k}^2}\,,$$
where  $w_k$ give less weight to less important frequencies. In sound, you can find A-weighting curves (and B-, C-, D- and Z-).

In images, a Matlab code function ratio = wsnr(orig, dith, nfreq) gives references to ancient papers:

J. Mannos and D. Sakrison, "The effects of a visual fidelity criterion on the encoding of images", IEEE Trans. Inf. Theory, July 1974.
T. Mitsa and K. Varkur, "Evaluation of contrast sensitivity functions for the formulation of quality measures incorporated in halftoning algorithms", IEEE Int. Conference on Acoustics, Speech and Signal Processing, 1993


How to find correlation between two signals regarding peaks? In a super-resolution task, I try to predict the higher frequencies of a signal.
X low-resolution (lr) is my input and Y high-resolution (hr) is the target. Where X contains only frequencies till 2 kHz and Y till 4 kHz.

I would like to know which frequencies in Y have a correlation with frequencies in X.
How can I find peaks in both areas (lr and hr) which occur always together?
for example, a peak at 1kHz always means there is also a peak at 3kHZ. In general, I would like to find patterns in the low-resolution input that helps me to find a correlation in high-resolution.

My low-resolution signal is generated from the high-resolution signal. It does not contain aliasing, because it was low-pass-filtered and downsampled correctly.

As an example plot : the lr and hr spectrum. Lr is generated by low-pass-filtering the hr spectrum. I would like to predict the higher frequencies.The lr image does not contain aliasing. So i need to find patterns in the lr signal for example maybe peaks, which always occur together with peaks in the hr signal. That would mean that there is a correlation between lr and hr.

Your premise is due to the physical nature of the device under test it is possible that occurrence of higher frequencies vibrations can be predicted by lower frequencies and this is important as it will allow you to use a lower bandwidth sensor that can only respond to lower frequencies.

Interesting problem and I would suggest this approach to find evidence of any such correlation:

You will be simply correlating “waveforms” to determine this and in this case the suggested waveform to correlate would be frequency versus time. You can establish these waveforms using a sliding windowed DFT. Decide on a window length (further details on that below) and starting at the beginning of your much longer data captures window a block of the data (more details below) and then use the FFT to compute the DFT for that block. Save the data as time index n=0 and then shift the window by one sample and repeat to get the data to time index n=1. Repeat this to the end of the data set such that the result for each DFT bin represents frequency versus time (a form of the STFT).  Once you have this result with frequency bins as columns and time index as rows you can cross correlate all columns to assess if there is anything significant. Anything that does not stand out significantly represents a noise sample, and you can rss all noise samples to get the noise floor of your measurement approach. I would say that anything that is <4 sigma is not significant and anything <2 sigma would reasonably be noise only. You can also do a histogram of the results which should be a standard bell curve but will also help show significant outliers.

Considerations for window size and window choice: with the length of the window you are trading frequency resolution with time resolution: the wider the window the narrower will be the frequency resolution. Specifically if you used a rectangular window the equivalent noise bandwidth of the result is 1 bin wide so the resolution in this case is $f_s/N$ where N is the number of samples or bins in the DFT and $f_s$ is the sampling rate. However the result will be the average frequencies over the time duration of the window so the trade space between time and frequency should be clear.

Next comes the window choice. With a rectangular window (no further scaling of the data prior to taking the DFT) the equivalent noise bandwidth is one bin as mentioned above. This is actually the best resolution you can achieve with the approach I describe as any other window will make the resolution wider. However most importantly is that with the rectangular window the shape of the equivalent filter in frequency is a Sinc function centered on each bin with first null at 1/T or 1 bin away. This means that there will be significant spectral leakage from other frequencies resulting in a false “correlation” if further windowing is not used. Windowing provides a trade of frequency resolution with dynamic range. The rectangular window has the best resolution (narrowest main lobe in its frequency response) but poorest dynamic range (highest sidelobes in its frequency response) compared to all other windows. It is important given your goal to window so that you don’t get a false sense of correlation due to those high sidelobes. I recommend the Kaiser window since you can adjust its Beta term to dynamically make the trade of resolution and dynamic range.

My guess is this would show that there is no significant correlation beyond integer harmonics but without having the mechanical model I can’t be certain of that so this will be interesting to see. If you do proceed with this,please post your results of what you find.


Question on LSI system Why second option (1) is not correct ?

Given an impulse response of
$$h(t) =
\begin{cases} \frac{1}{t^2} &,& \text{ for } t < 1/2 \\ 0 &,& \text{ otherwise } \\ \end{cases}
$$

Its associated CTFT (frequency response) is
$$H(\Omega) = \int_{-\infty}^{\infty} h(t) e^{-j \Omega t} dt $$

Now look at its value at $\Omega = 0$ which is $H(0)$;
$$H(0) = \int_{-\infty}^{\infty} h(t) dt = \int_{-\infty}^{\infty} \frac{1}{t^2} dt = 2 $$

Conclusion: since $H(0)$ is $2$ , then $1$ or $1/2$ or any number $x$ less than $2$; i.e., $x < 2$, cannot be an upper bound for $H(\Omega)$ as it already passes that value...


What happens when N increases in N-point DFT I am curious about DFT, and I wrote a simple MATLAB code to test what happens when $N$ increases. I took a rectangular signal with length $L=15$, an then found th DFT of 16, 32 and 64 points. I looked at the plots of the DFTs, however I could not understand why there is difference between them and how zero padding affects the DFT. Could someone explain it?

The length N of the DFT is the number of frequency points that will result in the DFT output. Zero padding will result in more frequency samples, however this does not increase frequency resolution, it just interpolates samples in the DTFT.   The frequency resolution is given by $1/T$ where T is the time length of your data (regardless of sampling rate). So if you want to increase the actual frequency resolution, you need to increase the number of samples at a given sampling rate, or decrease the sampling rate which would increase the time length for the number of samples you have.

For further explanation on the difference between the DTFT and DFT, see my post response here:

For 2D signals can it be said that the frequency response is the same as the Fourier transform?

Here is a demonstration showing frequency resolution and zero padding. The two red squares on the frequency axis at the top of the picture denote the "true" frequency location for two test tones, at f1 = 0.26 and f2 = 0.28.  Given the seperation of 0.02 in normalized frequency (cylces/sample) we would need to have greater than 1/0.02 = 50 samples to resolve the two frequencies.   The top plot shows the result of having only 10 samples of the data and zero padding out the time domain data for final datasets of 10 (no padding), 20, 50 and 100. We clearly see in this plot that as we add more samples, we interpolate more frequency points, but it does not offer any further information about the individual frequency content of our two tones.

Below that is another plot where instead of just adding zeros, we increase the number of samples in our dataset with the same two frequency tones present and again take the DFT with 10, 20, 50 and 100 samples.  At 50 samples we have enough frequency resolution for our separation in frequency (but still difficult to see on the plot), while at 100 samples the separation is clearly visible.

Why is this? Padding with zeros is the same as multiplying the longer sequence with a rectangular window.  Multiplying the time domain sequences is the same as convolution in the frequency domain  (circular convolution for the DFT).  The time domain rectangular window that we effectively used for the example above, and the resulting DFT of the window itself is shown in the plot below. The DFT of the window is what would convolve with our original frequencies. As the time length of the window increases, the width of the main lobe of its transform in frequency gets more narrow (the nulls are located at 1/T) and therefore the frequency resolution increases.


What is impulse response (or equivalently the frequency response) of the channel ? Please allow me to ask this silly question. I googled, but I still not clear about this. What is impulse response of a channel (wireless communication channel). Someone can explain me with a lot of pictures? I know this question is too ez, but I really need to understand it fully. Thank you again.

This is a fundamental topic and it needs a lot of explaination.
I suggest you use the fading channel object of MATLAB, especially the Channel Visualization Tool.

There are lots of examples with visualization in this link. MATLAB channel OBJECT


Effects of linear interpolation of a time series on its frequency spectrum Situation
In order to synchonisize different time series i have to apply linear interpolation on them. After the interpolation and synchronization the
signal is transferred into its frequency domain for further analysis.
The interpolation step should therefore not alter the frequency domain information.
Experiment
Via python I created a signal of white noise.
The second signal is based on the first but shifted by half of the original sampling interval. The third signal was created by oversampling the original signal by the factor 5.0.
Afterward all three signals were transferred to the frequency domain (Bottom of the Figure)
Questions

Why does the curve of the shifted signal decline for higher frequencies?
Why does the curve of the over sampled signal decline for higher frequencies? Though the over sampled signal has the exact same shape as the original signal.
What can i improve in order to interpolate the original signal at specific points but minimize the effects on its frequency spectrum?

My python code for the experiment and the plot: https://pastebin.com/iz08Hiud
Many thanks in advance!
Update
Thank you all for your comprehensive adn helpfull answers, especially robert bristow-johnson and howpow2 for directing me on the right path. I am using now an implementation of the sinc kernel interpolation, which shows an almost perfect frequency response.
See the bottom row of the plot for the results in form of the error between the original frequency spectrum and the frequency spectrums of the interpolated signals

Duane Wise and i wrote a paper back in the 90s that we presented to an AES convention that spelled out how to model time-domain polynomial interpolation (of which linear interpolation is an example) in the frequency domain.

i think you can get a copy here: Performance of Low-Order Polynomial Interpolators in the Presence of Oversampled Input


Downsampling and Upsampling Using FFTW I'm having 2 signals which have different sampling rates i.e., 1ms(A) and 4ms(B) and I've tried to upsample/downsample either of the signals based on the code snippet Resampling based on FFT which I believe is equivalent of what is descriped in How to resample audio using fft or dft

% FFTRESAMPLE Resample a real signal by the ratio p/q
function y = fftResample(x,p,q)

% --- take FFT of signal ---
f = fft(x);

% --- resize in the FFT domain ---
% add/remove the highest frequency components such that len(f) = len2
len1 = length(f);
len2 = round(len1*p/q);
lby2 = 1+len1/2;
if len2 < len1
  % remove some high frequency samples
  d = len1-len2;
  f = f(:);
  f(floor(lby2-(d-1)/2:lby2+(d-1)/2)) = [];
elseif len2 > len1
  % add some high frequency zero samples
  d = len2-len1;
  lby2 = floor(lby2);
  f = f(:);
  f = [f(1:lby2); zeros(d,1); f(lby2+1:end)];
end

% --- take the IFFT ---
% odd number of sample removal/addition may make the FFT of a real signal
% asymmetric and artificially introduce imaginary components - we take the
% real value of the IFFT to remove these, at the cost of not being able to
% reample complex signals
y = real(ifft(f));

i.e., reimplemented the fftResample in C++ using FFTW libraries.

The problem is that the output is distorted because of the zero-padding.
I'm relatively new to Singal processing and I would appreciate if users can help me understand:

if I'm using the right code snippet to upsample/downsample?
highlight the benefits of upsampling(A to 4ms) compared to
    downsampling(B to 1ms)?

When you downsample a signal you must make sure that the Nyquist condition remains satisfied.

For example, say that signal $x$ has frequencies between $0-25$ Hz, and that the sampling frequency of $x$ is $100$ Hz. The Nyquist frequency is $100/2 = 50$ Hz, thus you could downsample $x$ by a factor of two with no problems.

Now, say that you would like to downsample $x$ by a factor of $4$ instead of $2$. This would mean that the new sampling frequency would be $25$ Hz. Since the new sampling frequency is $25$ Hz, it must not contain frequencies above $25/2=12.5$ Hz. However, $x$ had frequencies between $0-25$ so there would be aliasing. In order to avoid that, you must low-pass filter $x$ first between $0-12.5$ Hz.

The process of low-pass filtering + downsampling is called decimation. You can find more information about it on Wikipedia.

Upsampling would be basicly interpolation, so I would avoid that where possible. I did not look at the URLs you provided, please phrase your problem in full here.


Running a sox spectrogram on a bandpass filtered wav file Let's apply a 400-2000Hz bandpass filter on a respiratory .wav file:

sox audio1.wav audio1.bandpass.wav sinc -t 10 400-2000

Now I'd like to generate a spectrogram that takes advantage of the reduced bandwidth to still generate a full size image. But how to to do that? Notice what the result is from running spectrogram on the bandpass file:

sox audio1.bandpass.wav -n spectrogram -r -o audio1.bandpass.png -m

That's clearly a big waste of image real estate. What step(s) am I missing here to use the full image size to focus on that already-limited frequency bands?

Update From accepted answer: works great

Add a step after the initial bandpass:

  sox audio1.bandpass.wav -r 4000  audio1.bandpass1.wav

Then

sox audio1.bandpass1.wav -n spectrogram -o audio1.bandpass.png -m

sox is doing exactly what you tell it to do – it plots a spectrogram of the bandwidth you can represent with your sampling rate. The fact that most of your picture shows no energy just shows your filter is working, and that's exactly what one would want to see!

However, as you can see, your signal is solidly oversampled: You can downsample it to a much lower sampling rate after filtering! (and downsampling by a factor $N$ here really is just throwing away $N-1$ samples, keeping one, throwing away $N-1$ samples..., since your filter is an excellent anti-aliasing filter.)


Implementation of goertzel algorithm in C I am implementing BFSK frequency hopping communication system on a DSP processor. It was suggested by some of the forum members to use Goertzel algorithm for the demodulation of frequency hopping at specific frequncies. I have tried implementing  the goertzel alorithm in C. the code is follows :

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
int k,i;
float   floatnumSamples;
float   omega,sine,cosine,coeff,q0,q1,q2,result,real,imag;

floatnumSamples = (float) numSamples;
k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
omega = (2.0 * M_PI * k) / floatnumSamples;
sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;
q0=0;
q1=0;
q2=0;

for(i=0; i<numSamples; i++)
{
    q0 = coeff * q1 - q2 + data[i];
    q2 = q1;
    q1 = q0;
}
real = (q1 - q2 * cosine);
imag = (q2 * sine);
result = sqrtf(real*real + imag*imag);
return result;
}

When I use the function to calculate the result at specific frequencies for a given dataset, I am not getting the correct results. However, if I use the same dataset and calculate the goertzel result using MATLAB, then I get the results perfectly. I am implemented the algorithm using C, with the help of some online tutorials that I found over the internet. I just want to get the view of you guys if the function is implmenting the goertzel algorithm correctly.

Thanks,
Anshu

It's been a while, but if anyone would like a working fixed-point implementation (no floats/doubles) in C you might take a look at: http://www.ti.com/ww/cn/uprogram/share/ppt/c6000/Chapter17.ppt
– Slide 14.

It looks like it will work beautifully and efficiently in even a low-powered micro-controller.
(In which case some of the ‘int’ variables in that TI C code may need to become ‘longs’ to be 32-bit.)


How are the constellation point coordinates determined in digital modulation I need to understand where the concept of the distance comes from. What I mean is, how do we know what coordinates are to assign to points in the constellation. For example here for QPSK the points are given coordinates $$\pm\sqrt(E_b),\pm\sqrt(E_b)$$
For BPSK on the other hand they are $$\pm\sqrt(E_b),0$$
This question was answered here as:
"A symbol with coordinates $(x,y)$ has energy $x^2+y^2$", and then it makes sense as for example for BPSK the energy per symbol is just $E_b$ and $\sqrt(E_b)^2+0^2=E_b$ but where does this come from? How to derive this? Thanks.

Constellation diagrams exist in what is called signal space which is an abstraction used to describe finite-energy signals. The coordinate axes, even if they are marked $x$ and $y$ as in Marcus Muller's answer, really represent unit-energy signals such as $$s_I(t) = \sqrt{\frac 2T}\cos(2\pi f_c t)\mathbf 1_{t \in [0,T)}(t),$$ $$s_Q(t) = -\sqrt{\frac 2T}\sin(2\pi f_c t)\mathbf 1_{t \in [0,T)}(t)$$ where $T = \frac{n}{f_c}$ and $\mathbf 1_{t \in [0,T)}(t)$ means a rectangular pulse of duration $T$. Note that both sinusoidal pulses consist of and integer number $n$ periods of a sinusoid of frequency $f_c$. Verify for yourself that $s_I(t)$ and $s_Q(t)$ are indeed unit energy signals and that they are orthogonal, that is, verify that
$$\int_0^T [s_I(t)]^2 \mathrm dt = \int_0^T [s_Q(t)]^2 \mathrm dt = 1;\quad
\int_0^T s_I(t)s_Q(t) \mathrm dt =0.$$
If we use $\pm As_I(t)$ as our PSK signals, then the signal energy is $A^2$(why?) and so if we use $E_b$ to denote the energy per bit, then
with respect to the signal space with axes $s_I(t)$ and $s_Q(t)$, the two signals used in PSK can be represented by the constellation points $(\sqrt{E_b},0)$ and $(-\sqrt{E_b},0)$. Left-handed folks might prefer to use PSK signals $\pm As_Q(t)$ in which case the signal constellation would have points $(0,\sqrt{E_b})$ and $(0,-\sqrt{E_b})$.
Mealy-mouthed ambidextrous people might want to use
$$\pm A[s_I(t)\cos(\theta)+s_Q(t)\sin(\theta)] = \pm A\sqrt{\frac 2T}\cos(2\pi f_c t+\theta)\mathbf 1_{t \in [0,T)}(t)$$
as PSK signals in which case in the signal space with axes $s_I(t)$ and $s_Q(t)$, the constellation points would be $\pm (\sqrt{E_b}\cos\theta,\sqrt{E_b}\sin\theta)$. But if we had chosen axes $\sqrt{\frac 2T}\cos(2\pi f_c t+\theta)$ and
$-\sqrt{\frac 2T}\sin(2\pi f_c t+\theta)$ instead (this just corresponds to a rotation of axes by $\theta$ from the previous signal space), then the constellation points would once again be $(\sqrt{E_b},0)$ and $(-\sqrt{E_b},0)$.

Thus, in signal space in which the coordinate axes represent unit-energy orthogonal signals, the energy of the signal corresponding to a constellation point equals the squared distance from the origin. Note that in (antipodal) PSK considered above, this squared distance from the origin is $E_b$ for each of the two constellation points. For QPSK (which I have ignored for the most part but you can read about it  in this answer of mine), each of the four constellation points are at squared distance $2E_b$ from the origin which makes sense in that $E_b$ is the energy per bit and a QPSK signal carries two bits.

Finally, with regard to the importance of distance as a concept, given two constellation points in signal space, the probability of error (in an AWGN channel with two-sided power spectral density $N_0/2$) is $Q\left(\frac{d}{\sqrt{2N_0}}\right)$ where $d$ is the distance between the two points. For the PSK examples considered above, $d = 2\sqrt{E_b}$ which gives the familiar result
$$P_e = Q\left(\frac{d}{\sqrt{2N_0}}\right) = Q\left(\frac{2\sqrt{E_b}}{\sqrt{2N_0}}\right) = Q\left(\sqrt{\frac{2E_b}{N_0}}\right)$$


Matlab FIlter coefficients to CMSIS api I am trying to implement iir Band Pass filter design in STM32F429i discovery kit.
I konw that it needs coefficients to be generated from matlab and feed them into cmsis api to achivi filtered data.

My SOS Matrix look like shown below:

1,  0,  -1,  1, -1.9974,   0.99743,
1,  0,  -1,  1, -1.9995,  0.99947,
1,  0,  -1,  1,  -1.9937,  0.99371,
1,  0,  -1,  1, -1.9982,  0.99817,
1,  0,  -1,  1, -1.995,  0.99499,

above matrix i put into my coeff array removing column 4th which contains all 1's since its A0 in coeff.

static float32_t filter_coeffs[] = {
1,  0,  -1,  -1.9974,   0.99743,
1,  0,  -1,  -1.9995,  0.99947,
1,  0,  -1,  -1.9937,  0.99371,
1,  0,  -1,  -1.9982,  0.99817,
1,  0,  -1,  -1.995,  0.99499,
};

here is the code snippet am using :

#define NUMSTAGES 5
static float32_t firStateF32[2*NUMSTAGES];
arm_biquad_cascade_df1_init_f32(&S, NUMSTAGES, filter_coeffs, firStateF32);

when I am using matlab my signal are filtered perfect, but in microcontroller output is initially oscillating and than "1.#QNAN"
am i missing any step??
please correct me.

Regards,
Lokesh Bhatt

I know this is an old question, but maybe this is useful to someone.
In CMSIS documentation, in Biquad Cascade IIR Filters Using Direct Form I Structure, says that the used equation is:

$$y[n] = b_0 x[n] + b_1 x[n-1] + b_2  x[n-2] - a_1 y[n-1] - a_2 y[n-2]$$

You need to look closely to coefficients $a_1$ and $a_2$. You need to change the sign.

So in your case, if matlab gave you this result (SOS Matrix)

1,  0,  -1,  1, -1.9974,   0.99743,
1,  0,  -1,  1, -1.9995,  0.99947,
1,  0,  -1,  1,  -1.9937,  0.99371,
1,  0,  -1,  1, -1.9982,  0.99817,
1,  0,  -1,  1, -1.995,  0.99499,

then you should change it to:

static float32_t filter_coeffs[] = {

1,  0,  -1,  1.9974, -0.99743,
1,  0,  -1,  1.9995, -0.99947,
1,  0,  -1,  1.9937, -0.99371,
1,  0,  -1,  1.9982, -0.99817,
1,  0,  -1,  1.995,  -0.99499,
};


Periodic signal checking when using $\Sigma$ Is the following signal periodic?

$$ \sum_{κ=-\infty}^{\infty}\left[\mathrm{rect}\left(\frac{t+2κ}{10}\right)\right]+\cos\left(\frac{π}{75}t\right)  $$

where rect is the rectangular signal

Let's define two new signals:

$$x_1(t)=\cos\left(\frac{π}{75}t\right)$$

$$x_2(t)=\sum_{κ=-\infty}^{\infty}\left[\mathrm{rect}\left(\frac{t+2κ}{10}\right)\right]$$

It's easy to see that the signal in your question can be calculated as

$$x(t)=x_1(t)+x_2(t)$$

The period of $x_1(t)$ is $T_1=150$. That can be seen easily: when $t=T_1$, the argument of the cosine becomes $2\pi$, so the function "resets".

Now let's see what happens with $x_2(t)$. Each term of the summation is a rectangular pulse of width $10$ centered at $t=-2\kappa$. As the magnitude of the shift factor (in this case, $2$) is less than the width, some pulses will overlap with others when performing the sum.

In fact, for almost every $t$, there will be exactly five rectangular pulses adding up. I say "almost" because there are special cases. When $t\in\mathbb{Z}$ and is odd, there will be only four pulses adding there, but there will also be two "half-pulses" corresponding to the value the $\mathrm{rect}()$ function takes at the discontinuity (every pulse will be discontinuos at two odd integers only), which is $0.5$. Because the values of $\kappa$ extend all over the natural numbers, it is valid to redefine $x_2(t)$ so that:

$$x_2(t)=5 \ \forall t$$

You can try that using Octave. I ran this script:

kMin = -500;
kMax = 500;
shiftFactor = 2;
tMax = shiftFactor*kMax + 10/2;
tMin = shiftFactor*kMin - 10/2;
t = tMin:tMax;
x2 = zeros(length(t),1);

for k = kMin:kMax;
  for ii = 1 : length(t)
    if t(ii) > (-shiftFactor*k - 10/2) && t(ii) < (-shiftFactor*k + 10/2)
      x2(ii) = x2(ii) + 1;
    elseif t(ii) == (-shiftFactor*k - 10/2) || t(ii) == (-shiftFactor*k + 10/2)
      x2(ii) = x2(ii) + 0.5;
    end
  end
end

plot(t,x2,'r','Linewidth',1.5);
xlim([tMin-1 tMax+1]);
ylim([0 5.5]);
grid minor;

The values kMin and kMax are the limit values $\kappa$ takes. Ideally, they should be infinite, but that's obviously not possible on a computer. shiftFactor is the constant that multiplies $\kappa$ in the summation. You can play a little bit with that number to see how the signal changes. I hardcoded the width of the pulses ($10$). I used the convention that the pulse equals $0.5$ in the discontinuity. If you plot the result, you get something like this:

Note that the ramps appearing by the sides do not happen in the infinite case. They show up here because the pulses that should be adding up for those values of $t$ were cut off when letting $\kappa$ take values between $500$ and $-500$ only.

So $x(t)$ is basically a cosine with an offset. Therefore, the sum of the functions we defined at the beginning is periodic. Even more, the period of $x(t)$ is the same as that of $x_1(t)$. Ergo,

$$T=T_1=150$$


LMS adaptive filter relatively delayed signal and reference inputs I'm using two microcontrollers to implement an adaptive LMS filter to filter out the noise from the signal. One is recording the noise and streaming the data to the other microcontroller which is recording the (noised) signal and using the noise reference, filtering it out.

How does the delay between those signal affect the filtering quality? Say the noise reference input is delayed for 100us/1ms/10ms/100ms/1s, how does that affect the filtering process?

The adaptive filter tries to emulate the assumed filtering process between the noise reference signal and the actual noise in the noisy signal.

If $n(t)$ is the actual noise in the noisy signal, and $n_r(t)$ is the noise reference, it is assumed that there's a linear filtering relationship between the two:

$$n(t)=(n_r*h)(t)$$

where $*$ denotes convolution (filtering), and $h(t)$ is some unknown impulse response, which should be approximated by the adaptive filter.

The practical effect of a delay of the noise reference depends on the impulse response $h(t)$. A delay of $\tau$ in the noise reference means that the adaptive filter needs to model $h(t+\tau)$ instead of $h(t)$. This can be a good thing if $h(t)$ is almost zero in the interval $0<t<\tau$, because then optimal use is made of the length of the adaptive filter. However, if there is significant energy in $h(t)$ in that interval, then, after delaying the noise reference, the adaptive filter must approximate a non-causal filter, which it will have a very hard time to do.

So in general, a delay of the reference signal results in a shifting of the impulse response that needs to be approximated by the adaptive filter. If this shifting is helpful or not depends on the amount of shift and the properties of the impulse response.


what is the suitable design Method to the filter? hi i'm going to design a low pass FIR filter for an EEG signal for the detection of eppileptic siezure and i want to know what is the suitable design Method to the filter? i need this FIR filter just to smooth the signal to use it for discret wavelet transform.
thanks

For a quick filter design without getting into trouble, I recommend avoiding any IIR structure and design an FIR using either least squares (firls in matlab) or Parks McClellan (firpm in matlab, remez in octave). There are plenty of answers here on how many taps you will need, (such as How many taps does an FIR filter need?) so you start with that for initial guidance and then refer to the specific help for firpm or firls for further design details using those tools.

Below is a slide I have on my typical design path for simple linear phase FIR filters. Starting with the specifications I estimate the number of taps, and then use that estimate to come up with (floating point) coefficients. I then confirm (using freqz in Matlab, Octave or Python typically) passband and stopband performance and adjust number of taps accordingly. Once the number of coefficients and values is settled, I then map to a fixed point implementation if needed using the guidance further given below and again confirm performance specifications. For final verification, especially if I have gone to fixed point, I do actual SNR testing with test waveforms through the filter at the maximum and minimum signal levels I would expect in my application using test waveforms. This is a good way to confirm that there are no disastrous overflow conditions or poor assumptions on quantization levels needed to maintain a target SNR directly at a component level test.

Design Tools:

Matlab: firpm, firls

Octave: remez, firls

Python: scipy.signal remez, firls

Online Tools:

Equiripple (Parks-McClellan): http://cnx.org/content/m12799/latest/

Least squares: http://cnx.org/content/m10577/latest/

Other tips:

Scale the output, let the filter grow the signal!  I think it was fred harris that had the view of analog filters attenuating signals out of band while digital signals grow the signal in band. I found that very insightful. To adjust the output level, scale after the filter; don't be tempted to adjust the output by scaling the input or scaling your coefficients (as it degrades SNR, output noise is accumulated noise from every tap and scaling coefficients increases their own individual error terms), beyond scaling to the proper precision setting as suggested below.  This is the reason for "extended precision accumulators"; they let the accumulated value from every tap grow and then once accumulated truncate the output in one location.

Coefficient bit width: A reasonable rule of thumb is to use at least 2 more bits for your coefficients than your datapath width if your datapath was sized based on SNR (or 2 more bits than your SNR requirement directly if your datapath is much wider than necessary).

You can also size the coefficient width based on your rejection requirement using 6 dB/bit (if you need a filter with 60 dB rejection, use 12 bits for your coefficients, with the coefficients scaled to fill that precision of course). Typically actual rejection achieved if limited by coefficient bit width is on the order of 5 to 6 dB / bit; I like to add at least the extra 2 bits for margin unless my design is severely resource constrained in that area.

Least Squares vs Parks McClellan (equiripple): In most of my design applications related to communications I have been more concerned with overall signal and noise than an absolute limit at any one given frequency, so typically use the Least Squares design method. That is the distinction between the two design choices (LS or P-McC) explained simply: Least Squares will minimize the error in a least squares sense, while Park McClellan minimizes the absolute error. Given the same number of taps, the total noise due to pass-band ripple and finite out of band rejection will be less with the Least Squares design approach when integrated over the entire waveform in time or over all frequencies, while the peak error at specific locations in time or frequency can be worst.


finding the output of these three systems when connected in series? I have just started to study Oppenheim's "Signals & Systems, Second Edition" and there is this easy-looking problem that has evaded me its solution for the past 48 hours:

consider systems $S1$, $S2$, and $S3$, where

$ x[n] \xrightarrow{S1} y_{(1)}[n] = \begin{cases}
x[n/2] & \text{if n is even;}\\
0 & \text{if n is odd.}\\
\end{cases}$

$ x[n] \xrightarrow{S2} y_{(2)}[n] = x[n] + (1/2)x[n-1] + (1/4)x[n-2] $

$ x[n] \xrightarrow{S3} y_{(3)}[n] = x[2n]$

To work out the result of $x[n] \xrightarrow{S1,S2,S3} y_{(1,2,3)}[n]$, I first tried to work out the result of $x[n] \xrightarrow{S1,S2} y_{(1,2)}[n]$ as follows:

$x[n] \xrightarrow{S1,S2} y_{(1,2)}[n] = \\
y_{(1)}[n] + (1/2)y_{(1)}[n-1] + (1/4)y_{(1)}[n-2] = \\ \begin{cases}
x[n/2] + (1/4)x[(n-2)/2] & \text{if n is even;}\\
(1/2)x[(n-1)/2] & \text {if n is odd.}\\
\end{cases}$

after that, feeding $y_{(1,2)}[n]$ to $S3$, we have:

$y_{(1,2)}[n] \xrightarrow{S3} y_{(1,2,3)}[n] = \\
y_{(1,2)}[2n] = \\ \begin{cases}
x[2n/2] + (1/4)x[(2n-2)/2] & \text{if 2n is even;}\\
(1/2)x[(2n-1)/2] & \text {if 2n is odd.}\\
\end{cases}$

which simplifies to $y_{(1,2,3)}[n] = x[n] + (1/4)x[n-1].$

but the solution manual says the output is in fact $x[n] \xrightarrow{S1,S2,S3} y_{(1,2,3)} = x[n] + (1/2)x[n-1] + (1/4)x[n-2].$ :(

What I want to know is, who is actually wrong: me or the manual? Thanks a lot in advance!

You are right and the manual is wrong.

Given $S_1, S_2, S_3$ and the respective input-output signals as below :

$$\\x[n] \rightarrow \boxed{S_1} \rightarrow v[n] \rightarrow \boxed{S_2} \rightarrow w[n] \rightarrow \boxed{S_3}\rightarrow y[n] \\$$

$\\$

$ x[n] \xrightarrow{S1} v[n] = \begin{cases}
x[n/2] & \text{if n is even;}\\
0 & \text{if n is odd.}\\
\end{cases}$

$ v[n] \xrightarrow{S2} w[n] = v[n] + (1/2)v[n-1] + (1/4)v[n-2] $

$ w[n] \xrightarrow{S3} y[n] = w[2n] \\ \\$

You would compute the output $y[n]$ in three steps as: (you have compined step 1 into 2 actually)

$v[n] \xrightarrow{S2} w[n] $

$
\begin{align}
w[n] &= v[n] + (1/2)v[n-1] + (1/4)v[n-2]  \\ \\
&=  \begin{cases} {
x[n/2] + (1/4)x[(n-2)/2] ~~~~~~ \text{if n is even;}\\
(1/2)x[(n-1)/2] ~~~~~~~~~~~~~~~~~~~~~~ \text {if n is odd.} }
\end{cases}
\end{align}
$

$\\$

and then, feeding $w[n]$ to $S_3$, we have:

$w[n] \xrightarrow{S3} y[n] = w[2n] \\$

$ y[n] =  \begin{cases}
x[2n/2] + (1/4)x[(2n-2)/2] & \text{if 2n is even;}\\
(1/2)x[(2n-1)/2] & \text {if 2n is odd.}\\
\end{cases}$

Hence (your) the solution is right:
$$y[n] = x[n] + (1/4)x[n-1]. $$

solution manual is wrong:
$$x[n] \xrightarrow{S1,S2,S3} y[n] = x[n] + (1/2)x[n-1] + (1/4)x[n-2].$$

You can verify this simply by putting some test signals, or even easier by writing a simple MATLAB/Octave/Python script to implement the $S_1,S_2,S_3$ system.


How exactly to implement the fractional part of a 'fractional equalizer'? I have written my own LSE equalizer (that works on the bit level), and literature I have read says that 'fractional' based ones are much better. So how exactly does one get 'fractional' bits? Is it as simple as resampling the softbit decisions to 2x, 4x or whatever of their bit rate and applying the equalizer like that?

Thanks

A typical ("non-fractional") equalizer operates at the symbol rate. That is, the inputs that you apply to the equalizer are generally soft decisions on what the unequalized received symbols are. A linear equalizer then takes this stream of soft symbol metrics and applies a linear filter that (hopefully) corrects for any non-ideal frequency response found in the path from modulator to demodulator.

The only difference with a fractionally-spaced equalizer (FSE) is the sample rate that it runs at: instead of the one sample per symbol that you see with a "standard" equalizer, an FSE operates at some multiple of the symbol rate. Typical values I've seen before are $T/2$- and $T/4$-spaced equalizers, which operate at twice and four times the symbol rate, respectively.

These samples should be straightforward to obtain; whatever processing steps in your system are used to generate the symbol decisions (e.g. when obtaining symbol timing synchronization) almost certainly operate at a minimum of 2 samples per symbol. Instead of decimating the receiver output to one sample per symbol, keep the higher sample rate intact through the equalizer in that case.


Discontinuity at the edge of Chebychev window I am using Chebychev window for its narrow main lobe. The problem of chebychev window is that it has discontinuities at the edge, and it seems that Taylor window solves this issue.

More detail:
http://de.mathworks.com/help/signal/ref/taylorwin.html
http://en.wikipedia.org/wiki/Window_function#Dolph.E2.80.93Chebyshev_window

I've searched around but I can't find any information on how to implement a Taylor window. Any information on taylor window or suggestions on fixing this issue of edge discontinuities would be very appreciated.

You've probably seen that at the bottom of the mathworks page on the Taylor window there are two references. But I guess it might be hard access these publications.

It's indeed difficult to find the formula online, but on the bottom of this page you can find a link to a (not yet released) Octave function implementing the Taylor window. I guess it's not too hard to distill the formula from the code. And here is a C implementation of the Matlab function taylorwin.m.


Removing noise from Signal captured by Laptop Sound Card using Mic In (4-pole Audio Jack) I am new to Signal Processing and I need a bit of help. I need to capture an Audio signal using the laptop Sound card. But before I captured the actual signal, I tested for the noise levels.

Generated a 1000Hz sine wave from a Mobile App and recorded the Signal using the laptop Sound card using a 4-pole audio cable on both sides (Mobile & Laptop).

After the Recording Audacity's spectrum analysis shows a lot of noise at all other frequencies.

A few questions:

Is there a way I can remove the noise levels from this captured signal since I know that this is a 1000Hz Sine wave?

My ultimate goal is to use the laptop Audio Card to capture the Signal and remove most of the noise.

Since I know the noise levels being introduced by the Sound Card, Can I Subtract a base line noise level from sound Card from the Signal to generate a "Clean" final signal?
Am I thinking about this in a wrong way being a newbee?

Thanks for your help.

If you know exactly which signal you want to record you could use a bandpass filter for this frequency.

This filter allows only a part of the signal to pass, like the name implies.
But otherwise, its hard to remove noise because this noise is over the whole spectrum.


How to get the channel impulse response by using PRBS signal with BPSK modulation? I want to get the channel impulse response by sending PRBS. Following is some channel measurement flowchart that I found

a(t) is the PRBS，sending it to transmitter to upconversion to BPSK modulated signal s(t)，r(t) represnet the receiving signal that pass throung real wireless channel. The receiver would downconversion the receving signal r(t) into I(t) and Q(t) signal.
After that, taking a(t) and I(t)+jQ(t) to do cross correlation to get the channel impulse response.

Yet, there is still some point that I confused about：

As far as I know, BPSK modulate and demodulate only need one I(t) signal，why the receiver output I(t) and Q(t) two signal?
After downconversion r(t), we get I(t)、Q(t) baseband signal. I am wondering aoubt the form of I(t)、Q(t) signal if I represnet a(t)、s(t)、r(t) as follow.

2

the point that I got stuck in is that is I(t) Q(t) signal BPSK demodulated signal or not? Moreover, no matter I(t) Q(t) is demoudulated signal or not, I both have problem. I describe it as below

Case 1： if I(t) Q(t) is BPSK demoudulated signal

The channel would definitely caused delay which would make the signal with phase shift that system unknown, so how can receiver demodulate the signal accurately?
Case 2：no need to do demodultion at Rx

To the best of my knowledge, if I have PRBS a(t) and the receving PRBS a(t-delay),then, after cross correlation of two signal I can get the delay time. Howerver, what if the receiver didn't demodulate the signal, how can I get the delay time and channel impulse response? Say the problem in another way is that, why calculate corss correlation of a PRBS signal and a sinusoid(I+jQ) signal can get the channel impulse response?

Do not think of your transmitter as BPSK since you are not intending to do data transmission (as you will not be doing data demodulation in the receiver) but simply an upconverter of a bipolar waveform (yes as far as the trasmitter operation goes, this is indeed identical to BPSK but I think that is giving you much confusion). So the transmitter is frequency shifting the baseband signal to an arbitrary carrier frequency.

The receiver simply down-converts the signal back to baseband.

Because of the phase and frequency offsets (in addition to the multipath distortion) the signal as received will no longer exist at baseband as a real signal but will have real and imaginary components. We need to maintain all of these components in order to solve for the channel. Therefore you must do a quadrature down-conversion that would maintain the real and imaginary components (I and Q) of the complex waveform. A very simple example is a channel with an impulse response represented by an impulse at t = 1 sec with a phase of 45°.

$$c(t) = \delta(t-e^{j\pi/4})$$

The received signal is the convolution of our transmitted signal with the channel, so in this case the BPSK transmitted signal would be received as a complex having a 45° phase rotation:

This is a very simple channel; in practice you will have multiple delayed copies at different delays and amplitudes combining in a possibly unrecognizable pattern in the receiver. If the transmitted waveform and the resulting received waveform is known, then you can use inverse convolution techniques such as the Wiener-Hopf equations to solve for the unknown channel, I describe this in detail in the following post:

Compensating Loudspeaker frequency response in an audio signal

(For the case of solving for the channel instead of the channel compensator, swap Tx and Rx as described in that post)

A PBRS (Pseudo-Random Binary Sequence) is an excellent choice for a training sequence as it has spectral content across the entire band. The equalizer can only determine its compensation at frequency locations where energy is present.


I/Q sampling with just one ADC Usually I/Q sampling is performed on two signals with a 90° phase shift using separate ADCs. Suppose I have only a single (fast) ADC available to perform I/Q sampling, which approaches do exist? Is this an entirely stupid idea?
A bit of background: I've been toying with the idea of building an amateur radio HF receiver. Yesterday the Raspberry Pi Pico microcontroller board was released. However I was disappointed to learn that it has only a single 12 bit, 500kS/s ADC, making it impossible to implement a conventional I/Q baseband receiver with dual ADCs.
There are two approaches I found:

Connect the in-phase and quadrature components to two different pins. Alternate between the components of the signal in the ADC. Then interpolate the signal to remove the phase shift between the samples.
Mix the signal with a higher intermediate frequency and filter it in such a way that successive samples are equivalent to the I, Q, -I, -Q components of the signal at the desired frequency. See the drawing and explanation on the first page of this paper:
"Ziomek, C., & Corredoura, P. (1995, May). Digital i/q demodulator. In Proceedings Particle Accelerator Conference (Vol. 4, pp. 2663-2665). IEEE." (PDF)

Is this an entirely stupid idea?

No, but you've just came to the conclusion that instead of sampling complex, with Nyquist rate being the bandwidth, you should do twice as many samples.
That simply means you're not doing IQ sampling, but low-IF or direct-RF sampling.

Mix the signal with a higher intermediate frequency and filter it in such a way that successive samples are equivalent to the I, Q, -I, -Q

so, multiply the real and imaginary parts with 1, -1, 1, -1 … each , then sample them.
1, -1, 1, -1 is is a cosine at exactly half the sampling rate. What you're mathematically doing is multiplying your IQ signal with a cosine + j· sine, each of half the sampling rate. $\cos ft + j\sin ft = e^{jft}$, i.e. that's a simple frequency shift.
That simply shifts up your signal from $f \in [-b/2; +b/2]$ to $f \in [0; b]$. I.e. you've just put your $b$-Bandwidth signal on an intermediate frequency of $b/2$. That's a low-IF heterodyne receiver, which existed before IQ was even discovered.
So, clearly, you can do that. I'd avoid shifting to only b/2, use something slightly higher if your bandwidth and sample rate allows, and forget all you've learned about complex sampling – we're back in real domain, and you simply have to sample with a rate that's higher than twice the bandwidth.

Also: As much hype there is around the RP2040: There's nothing new about having two cortex-m0 cores. It's not an especially fast microcontroller or anything, it really doesn't open any "new doors". The things that are nice about it are its flexible IO logic, but you need none of that.
In fact, these cores are ubiquitous, basically every semiconductor company sells microcontrollers based on that core. Although it's dual-core and clocked at a rather high rate, it's the last microcontroller that makes sense here: You need something that can process samples, not two high-clocked cortex-M0 cores that have no floating point and only rudimentary integer arithmetics. Get something with a Cortex-M4F instead. You might be getting a faster ADC, too.


Filtering using fourier series Suppose I have a measured signal $M$ that has frequency components from 0 to 50 Hz. I plot the specturm of this signal using FFT and I observe its frequency content (power vs frequency). Then, I decompose signal $M$ to its fourier series components and weights, I call this series $M_{hat}$. I compare the spectrum of $M_{hat}$ and $M$ and they look identical.

Now, I decide to extract the the frequency content between 0 to 10 Hz from $M_{hat}$. I consider the sins and cosines from 0 to 10Hz linearly combined with their weights. I call this signal $M_{hat10}$. When I plot the spectrum of $M_{hat10}$ and $M_{hat}$ on top of each other, I don't see them exactly matching in the range 0 to 10Hz. However, $M_{hat10}$ is indeed cut after 10Hz.

I don't understand why don't $M_{hat}$ and $M_{hat10}$ match in the range 0 to 10Hz.

This will match exactly if the time domain signal is properly recreated from the Inverse DFT as expressed in the sine and cosine form. In this case the DFT and inverse DFT are transform pairs so can be expressed in either domain with no loss of information.

Here is the general form for (painfully) computing the inverse DFT from sines and cosines which helps show the mathematical utility of always using the complex form of frequency when dealing with complex sinals ($e^{j\omega t}$ instead of $\cos(\omega t)+j\sin(\omega t)$). (And ultimately using the FFT algorithm and not the direct DFT equation is much simpler but this is insightful).

Given the Inverse DFT as:

$$x_n = \sum_{k=0}^{N-1}X_ke^{j2\pi nk/N}$$

When N is odd this is equally of the form:

$$x_n = X_0 + \sum_{k=1}^{(N-1)/2}\big(X_ke^{j2\pi nk/N}+ X_{[N-k]}e^{-j2\pi nk/N}\big)$$

If $x_n$ was real then $X_{[N-k]} = X^*_k$ leading to further simplifications, but in general for complex $x_n$ the above is given as below using $\omega_n = 2\pi n/N$ and the relationship $A e^{j\theta} = A\cos(\theta)+jA\sin(\theta)$:

$$x_n = X_0+\sum_{k=1}^{(N-1)/2}(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^{(N-1)/2}(X_k-X_{[N-k]})sin(k\omega_n)$$

When N is even the inverse DFT is of the form:

$$x_n = X_0+\sum_{k=1}^{N/2-1}\big(X_ke^{j2\pi nk/N}+ X_{[N-k]}e^{-j2\pi nk/N}\big)+ X_{[N/2]}e^{j\pi n}$$

Which similarly becomes:

$$x_n = X_0 +(-1)^nX_{[N/2]} \ldots \\+ \sum_{k=1}^{N/2-1}(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^{N/2-1}(X_k-X_{[N-k]})sin(k\omega_n)$$

In the OP's case where only the first $m$ frequencies of the DFT are used (corrresponding to $X_0$ when $m=0$, $X_1$ and $X_{[N-1]}$ when $m=1$ etc... and notably $X_{[N/2]}=0$, then regardless of $N$ is odd or even the result is:

$$x_n = X_0 + \sum_{k=1}^m(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^m(X_k-X_{[N-k]})sin(k\omega_n)$$

I also want to add this important additional note that this is actually a very poor way to implement a filter (it is called the frequency sampling approach, where the coefficients of the implemented FIR filter would be the IFFT of this frequency domain mask). It will indeed provide the results exactly as given in frequency at each bin center (in this case pass each frequency at the bin centers below 10 Hz with magnitude 1, and each frequency above 10 Hz with magnitude 0), but for all frequencies in between the bin centers both in the passband and stopband it will have a lot more deviation from other optimized filter design approaches (such as least squares and equi-ripple filter design).

This response at this post by @hotpaw2 explains this further in more detail:
 Why is it a bad idea to filter by zeroing out FFT bins?.


Simulate butterworth filter with initial condition of state vector in R I'm working in R with the packages signal and control (don't know wether these are even the best for my use case).
I have a high pass filter (butterworth, of order 1, some chosen cutoff frequency) and a (onedimensional) signal $s$.
I would like to filter $s$ using my high pass filter in a way, that pretends, that my high pass filter has previously seen the input $s[0]$ for an infinite amount of time.
I got a Python script, that does this using functions from the package signal. The procedure is as follows:

Create the filter (in tf-representation)
Switch to state space representation Calculate the initial condition for the state vector
Pass the filter, the input, the initial condition for the state vector to the function lsim2 of Python's signal package.
There you have the results

In R I can accomplish all but the second to last step.
Using the lsim function of control package doe not give the same results as in the Python scripts. I also tried the filter function from the package signal but this doesn't seem to allow me to provide an initial condition for the state space vectors.
Any hints on how to achieve this would be greatly appreciated.
I'm coming over from math.stackexchange so this is my first question.
Please be gentle, if due to my lack of experience here, the question is not perfectly posed.
Thanks in advance!

If you want to filter signals with an IIR filter, it's best to use the "standard" filter function of the language. That would be filter() for Matlab, lfilter() for Python and for R I believe this is signal::filter().
All of these functions allow you to specify the initial state. The exact details of how this is done, depends on the specific implementation of the filter.
If you don't specify any initial state the functions typically assume it's zero.
If I understand you correctly, you want to initialize the filter with the condition of $s[n] = s[0], n < 0$, i.e. a non-zero constant.
Apparently signal::filter() allows you to specify the arguments. initX and initY. I would assume the length of this arguments must match the filter order. You can set initX= s[0], for a first order filter. For an Nth order filter you would have to repeat s[0] N times. You would also have to calculate the steady state output for a constant input, but since it's a high pass filter, that should just be zero, so initY = 0 shuld work.
CAVEAT: I don't speak R and don't have an R development system to work with. This answer is derived from "best practices" and reading some documentation.


When is a LTI IIR $H(z)$ system minimum phase EDIT: I failed to mention that the system's inverse also needs to be causal and stable.
I cannot wrap my mind around on how when a system and its inverse are both causal and stable and LTI IIR  it is safe to say that the system is also minimum phase.
So the causality for a system means that  $$h[n]=0, \mbox{for all } n<0,$$ but what does the stability mean when talking about a system? That there are no poles inside the unit circle? and if so how does it turn out that when a system is causal and stable it's also min phase.
Also intuition-wise when is a system inverse causal and stable?
I should clarify that I'm talking about discrete time.

but what does the stability mean when talking about a system?

It means that all poles are INSIDE the unit circle.

and if so how does it turn out that when a system is causal and stable its also min phase.

Sorry, you got this wrong. Causal, stable and LTI does NOT imply minimum phase. A simple counter examples is a one-sample delay. It's causal, stable and LTI but definitely not minimum phase.
Minimum phase requires that all zeros are also INSIDE the unit circle as well.
Update
When you invert a system the poles become the zeros and the zeros become the poles. A minimum phase has both all zeros and poles inside the unit circle, that means that poles and zeros of the inverse will also be inside the unit circle and hence it's also minimum phase.
A non-minimum phase system has zeros outside the unit circle. If you invert the system, this zero turns into a pole that's outside the unit circle and hence the inverse is unstable (or non-causal).


temporal high-pass filtering to approximate phase derivative I am currently following the algorithm as based on a previously published research article "A Motion Detection Algorithm Using Local Phase Information" by Lazar et al. The basic overview of the algorithm is as follows:

I am able to do everything upto taking the 2D FFT. Once I take the FFT of a block of size 32 x 32 pixels, (Gaussian window has already been applied to this block), I get an array with the contents as follows:

For every entry in the 32 x 32 array, I understand that you can extract phase using the angle() function of MATLAB. However, I am not quite understanding

what the frequency(wx, wy) is referring to
how to temporally high pass filter phase information to approximate the change in local phase information.

Thank you.

I am not familiar with the algorithm you have proposed, but I believe I can answer the question you asked.

Once you have computed the DFT of your image, what you get is a 2D array of complex number. Each value in this array represent the amplitude and phase of a specific spatial-frequency (not time-frequency). Basically, what your paper tells you, is that you must loop every 1024 value into your 32x32 array.

Consider the output of the 2D DFT as a new image. You need to loop every pixel, extract the angle of the complex number and filter that angle over time, and you do this 1024 times.

I believe your algorithm applies to a video feed, if that is correct, each pixel instensity is located in a 3 dimensions reference system : $I(x,y,t)$. After the 2D DFT, you get new values (complex values) that are located in a new 3D reference system : $I'(\omega_x, \omega_y, t)$. If you look at a single point (or a "frequency pixel", if you allow me the term)  in your 2D arrray individually, you have a time-series of complex values. You could denote it $I'_{pixel}(t)$. Consider the angle of each single point like this $\theta(t) = angle(I'_{pixel}(t))$. This is the time series that you must filter over time.

At this point any discrete filter could do, you'll have to figure which one gives you best result. The filter will be something like :

$$ \theta_{filtered}(t) = b_1\theta(t) + b_2\theta(t-1)+...+a_1\theta_{filtered}(t-1)+...$$

Where you choose the $a_n$ and $b_n$ coefficients to get a performant high pass filter.

And the nice part about matlab is that you don't have to iterate over each point as it is designed to do matrix operation. You could go with something like the code below, assuming your software is running in real time and that you grab a new frame everytime you call the function

img = grab_new_frame()
img_dft = fftshift(fft(img)) % 2D spatial DFT
angle_map = angle(img_dft)
filtered_angle = angle_map *b1 + last_angle_map*b2 + last_filtered_angle*a1 % Where an and bn are scalars
last_angle_map  = angle_map
last_filtered_angle = filtered_angle

Good luck!


What will be the filtered output? I tried to solve this question from basic

Here is my work

Image 1

Image 2

But the correct answer is Option $(B)$.What is the mistake i am doing?

First note that:
$$ \cos(2\pi 50 t) \longleftrightarrow 0.5 \delta(f+50) + 0.5\delta(f-50) $$

$$\sin(2\pi 150 t) \longleftrightarrow 0.5 j \delta(f+150) -j 0.5\delta(f-150)$$

Hence the baseband spectrum is:
$$ X(f) = 0.5 j \delta(f+150) + 0.5 \delta(f+50) + 0.5 \delta(f-50) - 0.5 j\delta(f-150) $$

Then a shift to right by 100 Hz yields (without sampling weight):
$$ X(f-100) = 0.5\left[j \delta(f+50) +\delta(f-50) + \delta(f-150) - j\delta(f-250) \right]$$

Shift left by 100 Hz:
$$ X(f+100) = 0.5\left[ j\delta(f+250) +\delta(f+150) + \delta(f+50) - j\delta(f-50) \right]$$

Shift right by 200 Hz:
$$ X(f-200) = 0.5\left[ j \delta(f-50) +\delta(f-150) + \delta(f-250) - j\delta(f-350) \right]$$

Shift left by 200 Hz:
$$ X(f+200) = 0.5\left[ j\delta(f+350) +\delta(f+250) + \delta(f+150) - j\delta(f+50) \right]$$

The sum of all the shifts for $n=0,\pm 1,\pm2$ :
$$X_s(f) = 0.5\left[
j\delta(f+350) +\delta(f+250) + \delta(f+150) - j\delta(f+50) +
j\delta(f+250) +\delta(f+150) + \delta(f+50) - j\delta(f-50) +
j \delta(f+150) +\delta(f+50) + \delta(f-50) - j\delta(f-150)  +
j \delta(f+50) +\delta(f-50) + \delta(f-150) - j\delta(f-250) +
j \delta(f-50) +\delta(f-150) + \delta(f-250) - j\delta(f-350) \right]
$$

After the 100 Hz lowpass filter only those impulses inside $-100<f<100$ band remain:
$$LPF\{ X_s(f) \} = 0.5\left[ - j\delta(f+50) + \delta(f+50) - j\delta(f-50) +
\delta(f+50) + \delta(f-50) + j \delta(f+50) +\delta(f-50) + j \delta(f-50)
\right] $$

Which can be further simplied into:
$$LPF\{ X_s(f) \} =  \delta(f+50) + \delta(f-50)  $$

Hence we arrive (without sampling weight)
$$ \boxed{ z(t) = 2\cos(2 \pi 50 t)} $$

Therefore we can conclude that the output of the LPF will be proportional to $\cos(2\pi 50 t)$.


Will the received waves be the same if the transmitter use the beamforming to transmit the same signal? If there is a  transmitter with three transmitted antennas,and a receiver with five  antennas,and these five antennas are very close to each others.

Now transmitter use beamforming with these three transmitted antennas to transmit the same signal to the receiver,these five antennas will also receive the signal,but will their received signal wave be the same?i mean same magnitude,same phase.

Because if their phase are not the same,we should multiply different $e^{j\theta}$ to received signal that every antenna receive,then sum them,then decoding it.

But if their phase are totally the same,we can just sum them,them decoding it

Say you have the receive antennas very close together, like with spacing $\lambda/10$ or similar and they are uniformly spaced in a linear array. As you pointed out the the phase shifts between them will be fairly small, so what does it mean? It means that the spatial resolution of your beamformer will be poor. In the extreme case, consider something crazy where the spacing is $\lambda/1000$, all the receive antennas essentially receive the same copy of the signal. So why can't you just add them all up if they have no phase difference? You can, you just won't have any spatial discrimination and it is only really amplifying your signal by the number of receive elements.

I will leave a plot here that I generated in MATLAB which is the take away. You see as the spacing (specified in wavelengths) becomes smaller and smaller the beam response becomes flatter and flatter indicating that the receive array is unable to spatially discriminate because of the reasons outlined earlier. Hope this helps!


How Is the Gaussian Kernel Related to the Euclidean Distance of the Neighborhoods in the Non Local Means Denoising Algorithm? I am reading the paper A non-local algorithm for image denoising that describes the original non-local means algorithm. They define (p. 3 of the pdf) the distance between two square fixed neighbourhoods of two pixels as

\begin{align}
\|v(\mathcal{N}_i) - v(\mathcal{N}_j) \|_{2, a}^2
\end{align}

where $a> 0$ is the standard deviation of the Guassian kernel. However, I don't understand this notation. How is the Gaussian kernel used in the formulation? How is it related to the Euclidean distance?

I would appreciate if someone more familiar with image processing and computer vision and all this notation can give me an explanation.

So we have 2 neighborhoods (Which are Matrices) $ v \left( {\mathcal{N}}_{i} \right) $ and $ v \left( {\mathcal{N}}_{j} \right) $. How could we calculate the distance between them?
One option would be the Euclidean Distance:
$$ {\left\| v \left( {\mathcal{N}}_{i} \right) - v \left( {\mathcal{N}}_{j} \right) \right\|}_{2}^{2} = \sum_{k} {\left( {v \left( {\mathcal{N}}_{i} \right)}_{k} - {v \left( {\mathcal{N}}_{j} \right)}_{k} \right)}^{2} $$
Yet this gives each pixel in the neighborhood window the same weight. The writes of the Non Local Means Denoising Algorithm thought it would be better to give the pixels near the center of the window higher weight.
They did it easily using by multiplying the difference window by a Gaussian Function and then summing. Something like:
$$ {\left\| v \left( {\mathcal{N}}_{i} \right) - v \left( {\mathcal{N}}_{j} \right) \right\|}_{2, a}^{2} = \sum_{k} {w}_{k} \left( a \right) {\left( {v \left( {\mathcal{N}}_{i} \right)}_{k} - {v \left( {\mathcal{N}}_{j} \right)}_{k} \right)}^{2} $$
Where $ {w}_{k} \left( a \right) $ is a weight based on a 2D Gaussian function with its center aligned to the center of the neighborhood window and its Standard Deviation is given by $ a $.
The Weight Function
In the classic Non Local Means implementation the Gaussian functions is used as weighing.
Assuming the $ v \left( \cdot \right) $ operator is the Vectorization operator (Column Wise) then for a $ 31 \times 31 $ window the following MATLAB code will define $ {w}_{k} \left( \sigma \right) $:
windowRadius    = 15;
kernelStd       = 3.5;

vX = -windowRadius:windowRadius;
vK = exp(-((vX .^ 2) / (2 * kernelStd * kernelStd)));
vK = vK ./ sum(vK);
mW = vK.' * vK;
vW = mW(:);

figure();
surface(vX, vX, mW);

The weight function is given by:

Indeed centered pixels will have higher weight.