Design of 10th order high pass filter using 2nd order biquads with Matlab and further substitution on to TI AIC3204 I am involved in designing a 10th order high pass filter with cutoff frequency 300Hz.

My hardware AIC3204 has a 5 biquads of direct form I, 24bit coefficient
so I have chosen the FDA (filter design and anylysis) tool of MATLAB to design a 10th order filter and convert it to SOS sections 5 such and realize coefficients from it.

This is my procedure:

and realization of direct form I

Coefficients obtained as a file:

%
% Generated by MATLAB(R) 7.12 and the Signal Processing Toolbox 6.15.
%
% Generated on: 05-Aug-2013 15:26:23
%

% Coefficient Format: Hexadecimal

% Discrete-Time IIR Filter (real)
% -------------------------------
% Filter Structure : Direct-Form I, Second-Order Sections
% Number of Sections : 5
% Stable : Yes
% Linear Phase : No
% Arithmetic : fixed
% Numerator : s24,22 -> [-2 2)
% Denominator : s24,22 -> [-2 2)
% Scale Values : s24,23 -> [-1 1)
% Input : s16,15 -> [-1 1)
% Output : s16,11 -> [-16 16)
% Numerator State : s16,15 -> [-1 1)
% Denominator State : s16,15 -> [-1 1)
% Numerator Prod : s40,37 -> [-4 4)
% Denominator Prod : s40,37 -> [-4 4)
% Numerator Accum : s40,35 -> [-16 16)
% Denominator Accum : s40,35 -> [-16 16)
% Round Mode : convergent
% Overflow Mode : wrap
% Cast Before Sum : true

SOS matrix:
400000 000000 c00000 400000 7ab444 3b9018
400000 000000 c00000 400000 8bd4b8 377a90
400000 000000 c00000 400000 7404d3 34d7f8
400000 000000 c00000 400000 976d93 2b9955
400000 000000 c00000 400000 067eae d364fd

Scale Values:
77c3b6
77c3b6
6f5fe0
6f5fe0
6c9b03

Is my approach correct?

How can these higher order SOS filters be realised?

Your basic approach seems correct so far.
To implement your 2nd order filters you have to cascade the biquad filters, so that the output of the first feeds into the input of the next and so on.

This is due to the theorem that the overall transfer function $H(z)$ of $n$ cascaded subsystems is the multiplication of the transfer functions of the respective subsystems $H_i(z)$ like this

$$ H(z) = H_1(z) \cdot H_2(z) \cdot \,\, ... \, \cdot H_n(z)$$

So what Matlab gives you, are the coefficients $a_j$ and $b_j$ for the biquad filter transfer functions (called SOS matrix in your result):
$$ H_{biquad}(z) = \frac{b_0 + b_1z^{-1} + b_2 z^{-2}}{1 + a_1z^{-1} + a_2 z^{-2}}$$

To implement it you should use the Direct Form I in time domain - its just the inverse Z-transform of the transfer function above:

$$y[n] = b_0 x[n] + b_1 x[n-1] + b_2x[n-2] - a_1y[n-1]-a_2y[n-2]$$

Now you just have to cascade the filters, so that the output $y[n]$ of the first feeds into the input $x[n]$ of the next biquad and so on.

You may also want to check a similar question asked here: How does cascading biquad sections for higher order filters work?

And also the documentation of the Matlab biquad filter where it tells you how the SOSMatrix and ScaleValues are defined:
Matlab Biquad Documentation


Additive synthesis vs. "just draw the waveform" I've been experimenting with PYO recently, a Python module for sound synthesis and audio processing.  This is really my first foray into this topic, so I have what is probably a very basic question.

There are a few alternative mechanisms available for creating periodic waveforms such as a sawtooth or triangle wave:

Using additive synthesis (via SawTable, which internally calls HarmTable)

This approximates a waveform through successive addition of sine waves.  E.g.,

from pyo import *
table = SawTable(order=20)

Which gets me:

Using a "linear" table (LinTable)

This allows one to construct a waveform as a series of straight line segments.  E.g.,

from pyo import *
table = pyo.LinTable([(0, -1), (8191, 1)])

Which gets me:

These two waveforms sound similar, of course, but they are qualitatively different.  What is the technical difference between the two?  Is there a reason to prefer one mechanism over the other?

If your goal is to get a sawtooth wave, then additive synthesis is a bad way to go about it as a sawtooth wave is a superposition of an infinite number  of sine waves (in the continuous case).  Many instruments, particularly woodwinds, are well modeled by a relatively few sine waves.  The exact mixture determines the timbre and thus what instrument it sounds like, e.g. an oboe as compared to a clarinet.  Violins, however, generate a slip-stick motion that is more like a sawtooth wave and thus are harder to model accurately and efficiently with additive synthesis.

Unless you are going for an "electronic" sound, I would say using LinTable by itself would be a somewhat unusual thing.  However, using it as input to a filter would make more sense.  For example, if you did use a sawtooth wave to model a violin, you'd almost certainly pass it through a (among other things) low-pass filter to round it out a little and change the timbre.


What's a simple way to shift the pitch? Basically I want to perform an FFT to find the frequencies of the sound, then stretch the frequencies (for example a stretch by 2 means copy the value of the fourier coefficient X(n) to X(n/2) then use an inverse FFT. Will this work ?

To stretch the spectrum, you will also need to do some high quality (Sinc kernel) interpolation to fill in any bins that don't have an integer submultiple source.


Most common modern day Image Segmentation techniques I was doing some reading up on image segmentation techniques and I was wondering about modern day, state-of-the-art segmentation algorithms.

What current segmentation techniques that are a 'must-read', i.e. currently most commonly used in the community? Which techniques did you came in contact with and found most effective and useful (and for which application)?

I do not know about many segmentation techniques, but I have been dealing with structures that offer a "choice" of segmentation pieces that can be further examined to produce a satisfying segmentation.

Hopefully somebody else can write about some different state-of-the-art segmentation method that I don't know much about.

A small introduction as to why it is nice to have choices, or levels, of segmentation for a different image: segmentation is an ill defined problem. The ground truth does not exist: the desired results always depend on the user requirements and specifications. A nice quote:

  Even for a fixed image, there may be more than one "best" segmentation because the criteria defining the quality of a segmentation are application dependent. This motivated us to focus our research on image-partitioning techniques delivering "puzzle pieces" that can be used (...) to generate a segmentation satisfying specific user requirements

(P. Soille: Constrained Connectivity for Hierarchical Image Partitioning and Simplification (2008))

There are hierarchical structures, hierarchical image decompositions that propose image partitions with varying complexities. These structures are most simply represented as tree structures, where each node represents a region in the image. The idea with the structures:

the leaves of the structure are regions of a fine partition, or over-segmentation of an image (e.g. pixels, flat zones - connected regions of same intensity, watershed)
the links in the tree represent merging or union of (neighboring) regions and forming a more complex region, and are constructed in a way most likely to produce regions corresponding to objects (hopefully :) )
the complexity of the regions increases along each branch from nodes towards the root of the tree
every level of the tree (also, every cut) is a partitioning/segmentation of an image (coarser the closer to the root it is)
the root of the tree covers the whole image domain

The segmentation than consists of examining the regions and their proposed unions, to determine the regions in the tree or the cut of the tree corresponding to the required precision, or some known properties about object of interest, or other pre-defined user specifications.

Trees (that is, hierarchical image decompositions) with such properties are:

$\alpha$-trees
$(\omega)$-trees
Binary partition trees (aimed at segmentation, and object detection)

In addition to links to papers already provided, some more specific and less specific current papers on these kind of segmentation techniques:

Cousty, J., Najman, L.: Incremental Algorithm for Hierarchical Minimum Spanning Forests and Saliency of Watershed Cuts (2011)
Soille, P., Grazzini, J.: Constrained Connectivity and Transition Regions (2009)

(more practical):

Serna, A., Marcotegui, B.: Attribute Controlled Reconstruction and Adaptive Mathematical Morphology (2013)
Merciol, F., Lefèvre, S.: Fast image and video segmentation based on alpha-tree multiscale representation (2012)


Does Radio Communication Have to Account for the Doppler Effect? I know that radio communications take place in specific frequency bands, and that the doppler effect can cause shifts in frequency when there is relative motion between the source and destination objects. Do moving objects such as planes and rockets have to account for this, or is this movement negligible compared to the speed of light? How about communications with high velocity spacecraft and probes on other planets?

Does radio communication have to account for the doppler effect?

Yes.
Would be pretty terrible if that wasn't the case: RADAR wouldn't work!

Do moving objects such as planes and rockets have to account for this,

Yes. Phones in cars and trains, too. Your 5G NR phone of the future operating above 60 GHz, will have to do that at walking speeds, too.

or is this movement negligible compared to the speed of light?

You can calculate that yourself! With $f_0$ being the original frequency of your signal, and $c_0=3\cdot10^8\,\frac{\text{m}}{\text{s}}$ the speed of light, the Doppler shift as seen by an object travelling at speed $v$ is $\Delta f= f_0\cdot\frac{v}{c_0}$.
So, no, for most microwave communication (that's basically everything "broadband" that you do these days – TV, cellular, Wifi, drone video downlinks, the SpaceX video stream from a rocket), that's not negligible, when you move at high speeds.
And really, 200 km/h is pretty high speed, for example, for typical cellular comms (phones!). Say, your phone is using a band at 1.8 GHz, and you're sitting in a car going 100 km/h = 28 m/s. Then, your Doppler shift is 1.8·10⁹ Hz · 28/3 · 10⁻⁸ = 168 Hz. That doesn't sound so much, but it means that it takes but 1.5 ms until your phase has rotated by 90°, if my in-head calculations aren't off.
Since high-rate communications usually needs the phase to be consistent, you need to measure your channel way more often than that – so, more than 1000 times a second, your channel changes due to Doppler, if you don't correct the Doppler shift first!
And that, at a relatively benign 100 km/h (LTE does 300 km/h of bullet trains), at a relatively benign 1.8 GHz (modern WiFi works at 5.8 GHz).

How about communications with high velocity spacecraft and probes on other planets?

Of course. When you want to receive what was transmitted by a low-earth orbit satellite (the lower a satellite is, the faster it has to orbit around earth), you first have to correct for Doppler, before you can do much else.


Eigen Function of Linear Time Invariant (LTI) System Whenever I try to find the derivation of eigen signal of a LTI system, they always start with assumption that the input signal is complex exponential. Is there any other way to prove that this is the case starting from the basic equation of LTI system
$$\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau = \lambda x(t)$$
where $\lambda$ is the eigenvalue of the eigen signal

You can use the Convolution Theorem to show that $ \int_{-\infty}^{\infty} h \left( \tau \right) x \left( t - \tau \right) d \tau $ is equivalent to Element Wise multiplication in the Frequency Domain.

Hence for the case $ x \left( t \right) = A {e}^{s t} $ where $ s \in \mathbb{C} $ the output is just a multiplication by scalar of the input (The value of the System Transfer Function at the Delta point of the input).

See Eigen Functions of LTI Systems.


Sampling and Reconstructing Digital Signal Say that I have a digital signal (a bitstream actually) whose output changes every 100 mS.  I don't know exactly when it changes relative to my sampling, only that that the next bit to capture appears after 100 mS.  (I'm not capturing a square wave, but rather the signal might or might not assume a new value (either 1 or 0) after 100 mS.)  So I sample the signal every 25 mS.  How do I reproduce the bitstream from the sampled stream over some period of time...I can miss the first few bits of the capture, but after that I want to get every bit of the output?  I assume that there is a common algorithm for doing this.

By the way, the solution can all be done in software (Java / C++) and does not have to be real time...I can collect all samples and then post-process it to extract the signal.

Yes, there is. You have to study a broad field in digital communication systems which is called synchronization. Generally speaking, much of the synchronization principles are the same across different parameters. Specifically, you need to implement a timing synchronization (also called timing recovery or clock recovery) subsystem in a digital receiver. You can do it in either in open loop manner (feedforward) or closed loop method (feedback, PLL based). Once you start implementing it, then you will have specific questions for which you can return to SE time and again.

In this particular example, your clock frequency is all right (bits are changing every 100ms --- a constant difference from your sampling intervals). Ideally, the output of a matched filter is sampled at the end the sampling period. Since you don't know the ideal sampling instants (maximum eye opening), you have to use either the knowledge of the data (say, a preamble), or you can use characteristics of your bitstream in a non-data-aided manner.

Edit: If you are sampling at 25ms, that makes your samples/symbol equal to 4. If post processing is not an issue, then just use digital filter and square timing recovery for the whole preamble. the closed form expression gives you an exact estimate of timing phase offset.

Then, you will have to interpolate your data according to the estimate obtained above and downsample the signal from 4:1.


Commenting frequency spectrum of a sampled data the question below asked in my latest signal processing exam. I have no idea about it and i could not find any solution in literature. Could you please help me finding out the answer. Actually i believe that i am quite good at signal processing until this question.

A voltage signal is sampled digitally in the laboratory at three different sampling frequencies, and FFTs are performed on the sampled data. In all cases, 512 data points are taken. Here is a record of the observations:
‒ Data sampled at 115 Hz: A peak appears in the frequency spectrum at 2 Hz, with amplitude close to 3.0 V. Another peak appears in the frequency spectrum at 45 Hz, with amplitude close to 1.0 V.
‒ Data sampled at 280 Hz: A peak appears in the frequency spectrum at 48 Hz, with amplitude close to 3.0 V. Another peak appears in the frequency spectrum at 95 Hz, with amplitude close to 1.0 V.
-Data sampled at 500 Hz: A peak appears in the frequency spectrum at 185 Hz, with
amplitude close to 1.0 V. Another peak appears in the frequency spectrum at 232 Hz,
with amplitude close to 3.0 V.

Question: What are most likely the true frequencies (and their corresponding amplitudes) in the voltage signal.

Consider the 3V component first (note that aliassing only affects frequency, not amplitude, so the 3V component will have the same magnitude at whatever frequencies it is aliassed to):

For Fs = 115 Hz the 3V component shows up at 2 Hz, so it's real frequency could be any value which satisfies f = 115*k +/- 2, i.e.

  2 Hz
113 Hz
117 Hz
228 Hz
232 Hz
...

Do the same thing for Fs = 280 Hz, where we have f = 280*k +/- 48:

 48 Hz
232 Hz
328 Hz
512 Hz
608 Hz
...

And then again for Fs = 500 Hz:

232 Hz
268 Hz
732 Hz
...

We can see that the only frequency which corresponds in all 3 cases is 232 Hz, so it looks like we have a 3V component at f = 232 Hz.

You can now apply the same approach to identify the frequency of the 1V component.


Algorithm for detecting the time where the signal is above a threshold I need to detect the time window where a 1D-signal is above a certain threshold. If it dips below the threshold briefly I'd like to merge the two windows, if it dips below for a longer time, split them. The data are received digitally as array $S[t]$ with equally-spaced $t$ intervals and represent precipitation forecasts from a weather model.

I want to find $t_{on}$ / $t_{off}$ in signals like the following:

I know this isn't terribly hard and I can probably come up with an algorithm myself, but I thought this is a problem that smarter people than me have already thought about.

The algorithm will be implemented in software and run on general-purpose datacenter servers (no specialised hardware). It doesn't have to be crazy fast, "reasonable" is good enough (we have a few 100k time series to process at once every now and then; each series has probably less then 1000 points).

I'd welcome algorithm suggestions, Google keywords, book/article recommendations. I'm rather new to data processing, in case it isn't blatantly obvious.

Edit:
Thanks everyone for the great comments and answers. This question really couldn't have turned out much better. What I'm worrying about the most is how to handle pathological cases, like "oscillations" around the threshold.

My own idea was having separate on/off-thresholds, with $ thresh_{off} = f * thresh_{on} $ (with $f$ in the range of [0.5, 0.99]). I also really like Marten's idea of temporarily modifying the threshold value after the signal crosses the threshold in a hysteresis-like manner.

Are there other strategies to handle this issue that I should know about?

Are there other potential issues that I need to know about?

There is a book by Basseville and Nikiforov called "Detection of Abrupt Changes : Theory and Application" that they released to the public as a PDF several years ago (it's out of print, now, I believe).

That book looks at the basic CUSUM (cumulative sum) algorithm and how to choose appropriate thresholds for it.


Is it possible to recover QAM constellation of an undersampled OFDM signal? The OFDM signal I'm interested in is rated for a 204.8MHz sampling rate!

I'm limited to sampling at 20MHz to give me just enough resolution into the range of subcarriers I'm interested in, but I just want to see if I can recover the QAM constellation from the sub-carriers I'm interested in.

I understand OFDM very well, but I'm used to having the full sample rate!
Thanks for any insider comments.

Well, apply the definition of OFDM to a subcarrier range of an OFDM signal: it's still OFDM, in all letters of that acronym.

So, on a superficial level, yes, this works easily¹. Looking closer:

OFDM signals aren't usually rated for some sampling rate, the systems transmitting them just work at a specific rate, which mathematically coincides with the product of subcarrier spacing and numbers of subcarriers.

The good thing about OFDM is you don't have to worry about what the spectrum looks like in between the subcarrier "frequency points"; the orthogonality guarantees you've got no signal from other subcarriers on any subcarrier frequency, and the way perfectly synchronized OFDM works is by the DFT – which essentially samples the frequency domain.

Now, have you ever wondered why in e.g. a waterfall plot or on a spectrum analyzer, the spectrum in the unoccupied carriers at the edge doesn't actually decrease as fast as the (ideliazed) 6dB / decade of sidelobe power suggest?

We often forget that yes, as you say, OFDM systems are actually discrete systems based on the DFT; often, the DACs involved are actually internally oversampling, even. This implies a "folding" imaging scenario, where sidelobe power that would fall out of the (nominal) DAC's bandwidth will be folded back – not a problem, because it just adds power from other sidelobes to sidelobes; at the transmitter, the orthogonality still stands!

Now, at the receiver, we're happy that the edge carriers are unoccupied – because the moment where we're not sampling-clock-synchronized and frequency-synchronized yet, we'll have a lot of fun with circularly shifted-in sidelobes overlaying the tentative correct subcarrier frequencies at both ends of the spectrum (assuming we critically sample).

Your subset-of-subcarrier-restricted OFDM receiver runs exactly into that hardship: No "leeway" for frequency offset and clock error to work out. Note that in this case, the frequency offset isn't "worse" to deal with than in the "normal" OFDM receiver case. It's just that you're guaranteed to have aliases from the adjacent carriers and their sidelobes also ending up in your signal – and whilst these have the same PSD as in the full OFDM signal, your processing gains on the DFT structure of the signal during synchronizations are but a fraction of what the full system has!

So, as Jason R commented:

Make sure the analog side of your receiver has adequate filtering. In theory, the sidelobes of adjacent channels don't matter, in practice, your sampling clock offset makes them hurtful. At any rate, you don't want other carriers to alias into your signal.
Make sure your receiver has plenty of dynamic range headroom. Really - OFDM has notoriously nasty PAPR, and you don't want dynamically appearing intermodulation products (hint: in many hardware SDR receivers, sampling rate and RF frequency generation are somehow coupled, so sampling rate-rationally-related intermodulation products often appear – especially bad here, because these can map "unseen" carriers to your carriers.
Use additional digital low-pass filtering to enforce unoccupied carriers at the edge of your observed band; yes, that'll take energy away from the carriers you want to see (remember, sidelobes in OFDM do carry information energy!), but it'll make your synchronization life easier.
Adapt the synchronization algorithm to your carrier subset. This might or might not be easy or even possible, considering we know nothing about how your system does synchronization.
Equalization might become a lot harder, if the original system was designed to interpolate channel state between neighboring carriers (and at bandwiths> 200 MHz, we can be pretty sure it does so to minimize the need for pilots), since you've got only a very partial view on things.

Astonishingly, the QAM aspect together with the last point might be a show-stopper here: QAM requires relatively good equalization, and since your part of the OFDM signal might give you less than an equivalent part of the pilot information in the overall system, you might simply be channel-coherency limited (either not enough pilots in frequency direction, or not enough pilots over time). So, make a channel model first, figure out how often you'll see pilots, and estimate whether things can work first, before spending too much time on this.

¹ assuming you can still frequency-synchronize. Since you said you understand OFDM very well, I assume you've got practical knowledge of OFDM synchronization methods – but in short: most do work on a restricted bandwidth, with slight modifications; you need to assume 0 subcarrier coarse frequency offset, anyways. Note that I'd still assume that if implementing an OFDM receiver, the synchronization would be the hardest part to get right. Timing synchronization is another challenge, but gets even more system-specific; if your system uses Schmidl&Cox to detect both freq and time, then you're more or less lucky and just need to adjust tolerances (and sample counts) to make this work.


Single Pixel Camera - Compressive Sensing I work with a setup that measures high frequency fluctuations in light using a photodiode. We  steer the light over a sample as we measure these fluctuations.

I am familiar with compressive sensing where we only sample a few points to reconstruct the full signal. The sample is of course done in the same dimensions as the signal - for audio, the sampling is across time for instance. For an image, we randomly subsample in both x and y coordinates.

I came across a single pixel camera which employs a Digital Micromirror Device to project a pattern onto a scene which is then imaged with, again, a photodiode. Most of the descriptions I've read say that the each pattern would in effect sum the light from the scene and repeated scanning with random patterns can yield images with similar quality as raster scanning.

How are these images reconstructed? If we sum the light, we lose where the individual sums of that light came from. How does repeated acqusition solve this?

First I explain how compressive sensing is leveraged into imaging reconstruction  and then a little bit on how CS is deployed in an imaging hardware.
Compressive Sensing
For the sake of simplicity lets assume our image is 1D (i.e. a row of pixels). And assume the image you want to ultimately construct is $X_N$ with $N$ being length of it. The simple choice of course is to stack up $N$ sensors very close to each other and basically build up a 1D camera. However if you have only 1 image sensor, you need to either move the sensor in $N$ positions as if there are $N$ sensors or as you've done, steer light $N$ times and capture the pixel value. Now, the using compressive sensing you can do $M$ measurements ($M << N$) and still construct the image.
Now lets assume each pixel from pixel number 1 to pixel number $N$ as unknowns of the following system of equations. The right side, $Y$ is simple an vector that each element of it is sum of a number of pixels 1 to $N$, e.g.
$$Y_1 = c_{1,1}\times P_1 + c_{1,2}\times P_2+..c_{1,i} \times P_i+...c_{1,N}\times P_N$$
$$Y_2 = c_{2,1}\times P_1 + c_{2,2}\times P_2+..c_{2,i} \times P_i+...c_{2,N}\times P_N$$
$$ ...$$
$$Y_M = c_{M,1}\times P_1 + c_{M,2}\times P_2+..c_{M,i} \times P_i+...c_{M,N}\times P_N$$
$P$ values are pixels and $c$ values are some coefficients (actually those masks in single pixel camera). So if $N$ is 1000, the size,  $M$, of "measured" vector (as in CS literature), would be something around 100, that makes us, end up with a system of equation that is called under-determined. This constitute following, if we consider all $c$ as a matrix $C$ multiplied in signal $X$:
$$Y_M = C_{M,N}\times X_N$$
Now, if $M>=N$, the things will be very easy, you have a system of equations with $N$ unknowns and $M$ equations, easily solved (just multiply both sides with inverse or pseudo-inverse of $C$).
Now the power of compressive sensing is that it reconstructs signal even in case $M << N$. The assumption however is the signal is sparse, i.e. 99% of elements of $X$ are zero (or very close to zero). This is not true if you look at an image, unless you are in dark room rarely you see any pixel of an image to be zero, however there is another trick. If you take Fourier or wavelet transform of image, you will be supersize how many elements are very small (or zero) compare to just a handful of very large elements. Hence, we say images are almost sparse in transform domain. Now, if $X$ is our image, $\alpha$ will be our transformed image as following and $\alpha$ will be really sparse:
$$\alpha = \psi \times X$$ so obviosly, $X$ would be equal to $\psi^{-1}\times \alpha$. OK, if we substitute for $X$ we will have:
$$Y_M = C_{M,N}\times \psi^{-1}\times \alpha$$
Single Pixel Camera
Now, we have system of equations with $N$ unknowns but many of those $N$ unknowns are zero. Therefore we can use compressive sensing to reconstruct $\alpha$ and consequently $X$.
Now, how to multiply the pixel values with $C$? The way single pixel camera deals is they use micro mirrors. In our 1D example we need $N$ micro-mirrors, which can be selected to reflect light or not as if we are multiplying the light rays with either "1" (for reflection) or "0" (for no reflection). Reflecting and focusing all rays back on our single sensor is the addition operation. Each time we set some mirrors to reflect back and focus on a single sensor and measure the light, we take a measurement, i.e. an $Y_i$. After the "Measurement" which is mentioned sum of production, is digitalized by our sensor read-out we can carry out rest in the computer and run compressive sensing recovery, etc, but the first part all are analog.
The benefit of course is far fewer measurements $M<<N$ which naturally leads to more good stuff.
Now, how compressive sensing recovery algorithms solve these type of underdetermined sparse system of equations is another story. But this is pretty much how single pixel camera works.
Simulation
Simply make the following process in MATLAB or Python:

Make $M$ number of binary 2D arrays with size of your image, these will be your masks

Calculated and sum up element-wise multiplication of each mask (ith mask) with image and call it $y_i$

Stack up the masks into a 3D array and also the $y$s into a 1D vector

Reconstruct using CS (for this you need to use many available CS recovery algorithms)

There is an example code here which might help (Disclaimer: I have not checked if it works)
https://github.com/kurokuman/single-pixel-camera-simulation


What is the effect of reducing NFFT less than the signal length I am analyzing a signal in MATLAB by taking its fft. The singal has $200$ samples and I am currently taking nfft = 256. I understand that this would result in zero padding. Now, I want to try to reduce the nfft to 128, 64 etc. I wanted to ask what will happen to the fft.

Matlab says for Y = fft(X,n), If X is a vector and the length of X is greater than n, then X is truncated to length n. What would this actually result in. Can anyone explain, please?

Reducing the NFFT samples in your DFT will result is less resolution in the spectrum you'll be observing. The frequency resolution is defined as
$$
f_k = \frac{F_s}{N}k\quad \text{with}\quad k=0, 1, \ldots, N-1\tag{1}
$$
In your $N$-point Discrete Fourier Transform (DFT), decreasing $N$ increases the spacing between the $f_k$. Or put in another way, for a fixed sampling frequency $F_s$,

\begin{align}
N_2 &< N_1 \implies \Delta f_{k_2} > \Delta f_{k_1}\\
\text{Where}\quad f_{k_1} &=\frac{F_s}{N_1}k_1\quad\text{and}\quad f_{k_2} =\frac{F_s}{N_2}k_2\\
\implies \Delta f_{k_1} &=\frac{F_s}{N_1}\quad\text{and}\quad \Delta f_{k_2} =\frac{F_s}{N_2}
\end{align}

So say, you have some frequency components (i.e. peaks) you're observing in your spectrum with a high $N$, gradually reducing this value, you'll see that the peaks are getting less and less localized visually. Note that it is only visually, your DFT got your spectrum, you play with $N$ (i.e. NFFT) to visualise more out of your spectrum. Together with $F_s$, it determines the frequency resolution as defined in $(1)$.

EDIT:

To compute the energy, you sum the squared amplitudes spectrum values. Look for instance at the MATLAB code below for  the sine wave
$$
s(t) = 3\sin(2\pi 150t)
$$

Fs = 500;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 2000;             % Length of signal
t = (0:L-1)*T;        % Time vector
f = Fs*(0:(L/2))/L;

S = 3*sin(2*pi*150*t);

NFFT = 2^nextpow2(L);
NFFT1 = L;
NFFT2 = NFFT/2;
NFFT3 = NFFT/4;

Y = fft(S, NFFT1)/NFFT1;
Y2 = fft(S, NFFT2)/NFFT2;
Y3 = fft(S, NFFT3)/NFFT3;

EY = sum(Y.*conj(Y));
EY2 = sum(Y2.*conj(Y2));
EY3 = sum(Y3.*conj(Y3));

Your FFT have to be normalized and the energy computed as sum of the squared amplitude spectrum. For the sample code above you get:

EY = 4.5
EY2 = 4.49644474953643
EY3 =4.49832143842713

Which is from the amplitude of the sine $3$, giving ($3^2/2$). Also showing the effect of taking different $N$ values for the FFT on resulting energy.


Why is Frequency Shift Keying more immune to noise than ASK? Is it because noise manifests itself in the transmitted signal's amplitude?

Looking through some literature shows that a side by side comparison of FSK with ASK does not give FSK a performance advantage over ASK in the presence of noise (1); however you can do some nonlinear processing with FSK in the demodulator that does produce a signal that is more tolerant to noise than ASK.

In some instances of FSK modulation/demodulation schemes, you can force the receive signal to clip (i.e. put the signal through a Schmidt trigger) so the amplitude noise is essentially eliminated.  You can't do this with the ASK signal because this eliminates the information in the modulated signal.  For the FSK signal, the clipped signal still carries frequency information with a slight phase error due to the noise signal.  There is a threshold at which noise can cause enough phase error to affect the demodulation even after clipping the signal, but below that threshold, the signal is largely immune to noise.


Is Ideal LPF BIBO unstable? In one of other discussions : How to find frequency response, stability, and causality of a linear system?

I found a comment which was quite strong and definitely caught my attention.

  An ideal low-pass filter is an example of a system that is not BIBO
  stable even though its frequency response is bounded for all $f$

I am following the definition of stability as per here in wiki http://en.wikipedia.org/wiki/BIBO_stability

Can anyone give me a proof that ideal LPF can indeed be BIBO unstable?

Of course, ideal LPF with infinite gain can produce unbounded output. The question is restricted to LPF when gain is finite.

This answer is a response to a comment by the OP on on yoda's answer.

Suppose that $h(t)$, the impulse response of a continuous-time
linear time-invariant system, has the property that
$$\int_{-\infty}^{\infty} |h(t)| \mathrm dt = M$$ for
some finite number $M$.  Then, for each and every
bounded input $x(t)$, the output $y(t)$ is bounded also.
If $|x(t)| \leq \hat{M}$ for all $t$ where $\hat{M}$
is some finite number, then $|y(t)| \leq \hat{M}M$ for all $t$
where $\hat{M}M$ is also a finite number.
The proof is straightforward.
$$\begin{align*}
|y(t)| &= \left |\int_{-\infty}^\infty h(\tau)x(t - \tau)\mathrm d\tau\right |\\
&\leq \int_{-\infty}^\infty |h(\tau)x(t - \tau)|\mathrm d\tau\\
&\leq \int_{-\infty}^\infty |h(\tau)|\cdot|x(t - \tau)|\mathrm d\tau\\
&\leq \hat{M}\int_{-\infty}^\infty |h(\tau)|\mathrm d\tau\\
&= \hat{M}M.
\end{align*}$$
In other words, $y(t)$ is bounded whenever $x(t)$ is bounded.

  Thus, the condition
  $\displaystyle\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$
  is sufficient for BIBO-stability.

  The condition  $\displaystyle\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$
  is also necessary for BIBO-stability.

Assume that every bounded input
produces a bounded output.  Now consider the input
$x(t) = \text{sgn}(h(-t)) ~\forall~ t$.  This is clearly bounded,
($|x(t)| \leq 1$ for all $t$), and at $t=0$, it produces output
$$\begin{align*}
y(0) &= \int_{-\infty}^\infty h(0-\tau)x(-\tau)\mathrm d\tau\\
&= \int_{-\infty}^\infty h(-\tau)\text{sgn}(h(-\tau))\mathrm d\tau
&= \int_{-\infty}^\infty |h(-\tau)|\mathrm d\tau\\
&= \int_{-\infty}^\infty |h(t)|\mathrm dt.
\end{align*}$$
Our assumption that the system is BIBO stable means that $y(0)$ is
necessarily finite, that is,
$$\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$$

The proof for discrete-time systems is similar with the obvious
change that all the integrals are replaced by sums.

Ideal LPFs are not BIBO-stable
systems because the impulse response is not absolutely integrable,
as stated in the answer by yoda.  But his answer does not really answer the question

  Can anyone give me a proof that ideal LPF can indeed be BIBO unstable?

A specific example of a bounded input signal that produces an unbounded output
from an ideal LPF (and thus proves that the system is not BIBO-stable)
can be constructed as outlined above (see also my comment on the main question).


Smoothing a discrete function with point wise variances In Stark's Introduction to numerical methods [1970] it suggests fitting local polynomials to smooth discrete data before the Fourier transform using least squares minimisation. This gives the following rules for degree three five-point smoothing:

for the first two entries:
$$
g_1 = \frac{1}{70} (69 f_1 + 4 f_2 - 6 f_3 + 4 f_4 - f_5) \\
g_2 = \frac{1}{35} (2 f_1 + 27 f_2 + 12 f_3 - 8 f_4 + 2 f_5) \\
$$

for the last two entries:

$$
g_{N-1} = \frac{1}{35} (2 f_{N-4} - 8 f_{N-3} + 12 f_{N-2} + 27 f_{N-1} + 2 f_{N}) \\
g_{N} = \frac{1}{70} (- f_{N-4} + 4 f_{N-3} - 6 f_{N-2} + 4 f_{N-1} + 69 f_{N})
$$

and for all other entries:

$$
g_{n} = \frac{1}{35} (-3 f_{n-2} + 12 f_{n-1} + 17 f_{n} + 12 f_{n+1} - 3 f_{n+2})
$$

Would a weighted least squares approach following this idea be valid (what would the formulae be?) or are there more widely recognised approaches for smoothing functions with known variances?

Other then the endpoints, this is really just an FIR filter. You can easily assess the impact of smoothing on the spectrum by looking at the Fourier Transform of the impulse response of the filter . In this case it's a notch filter which seems like an odd choice.

In general smoothing is in many cases similar or equivalent to a lowpass filter operation and will in general decrease high frequency content of the signal. The "best" choice of low-pass will always depend on your specific application and requirements.


MR Segmentation to draw marker over the MR Dicom Images I am doing R&D over the MR Segmentation to draw marker over the Heart MR Dicom Images. My Question is that how I can draw a such marker over the dicom image that show in attach/link image. that Region of Interesr ROI read from the tag or its part of Pixel Processing.

Typically research papers handling DICOM standard images would probably extract the ROI using pixel processing methods. However, a lot of medical imaging companies like to implement their own DICOM tags(public/private) to keep track of such information for various purposes.

For example, General Electric keeps the following tag enabled on most MR images:

VerticesOfTheRegion found at (0028,9503), which is according to the DICOM standard, a private tag(odd location tags are private)

In short: Its implausible to tell you,how was the ROI extracted without having access to the DICOM metadata of the image in question.


What Is the Algorithm Behind Photoshop's "Black and White" Adjustment Layer? Can someone explain the algorithm behind the "Black and White" adjustment layer in photoshop?

I have to reproduce it using C++ for an application that emphasizes non-red/magenta(ish) pixels from an image (with a percent-like configurable tolerance), and this resource showed the behavior I'm expecting.

Still couldn't reproduce it, but I found a clue:

Each pixel is defined by up to two controls, one additive (RGB) and one subtractive (CMY).

Each (color) image is composed of RGB components. when you add (or reduce) a constant value to all pixels only in RED components you will see the effect equivalent to moving the RED tab towards the right, and same way reducing the RED component by a constant will have the reverse effect.

Like wise you can increment/decrement each component by a fixed value as described. If you increase/decrease ALL RGB components by a same value this will be equivalent to change in brightness (basically you are adding/removing WHITE color).

The Cyan, Blue, Magenta - corresponds to such transformation in CMYK color space. (But i guess, Blue in this color space corresponds to mix of Cyan and Yellow. So that's a bit tricky. The transformation is essentially same for all.

The last element Tint: { Hue and Saturation } corresponds to same operations but here, the images is first transformed in HSV model and then HUE and Saturation is added/subtracted independently.

I don't know about the exact relationships of markers of the dial to the corresponding numbers but that you can figure out by trying out practical values.


Having Glitches trying to use a Sinc Function Interpolation on a Signal, by chunks of 1024 samples I am trying to interpolate a signal using a third party lib that performs Sinc Function Interpolation on a signal. It works well if I input all the samples at once but introduce phase distortion if I try to perform by block of 1024 samples.

Some context:
I am trying to increase or reduce the playback speed of an audio loop.
The audio loop is a stereo interleaved samples with a total of 352,800 samples (176,400 stero frames).

It works flawlessly when I input all the 352,800 samples at once, but this can't be performed real time.
Instead, I try to feed it with chunks of 1024 samples as requested by the audio interface.

To do that, I first calculate how much I have to feed the interpolation function to get an interpolated result that fits this requirement:

// how much samples we need to prune in the original samples
let to_take = asked_size * playback_rate;

Then I pass theses samples to the interp function that spits out excaclty 'asked_size' samples, theses samples have been interpolated.

The interpolation function takes the following parametters:

samples: the samples buffer to be interpolated (1024*playback_rate)
sinc function (this is the sinc function as in the lib).
playback_rate
delay (which I set to 0.0)
wing (how many samples to looks from each sides for interpolation / convolution).

When the playback rate is 1, everything is smooth, but when I move this value, the global signal get distorded, because there is gaps between the first sample of the next buffer and the last sample of the previous buffer that I send to audio.

I suspect there is something to do with the phase, or something to do with the fact that the interp function miss some contexts on previous chunks of samples (blocks of 1024) for the outpout (concatenation of all the blocks of 1024) to be seemsless, and this results in nasty distorted audio.

How should I correct this behaviour ?

EDIT: Added implementation in Rust:

// block loop
pub fn next_block(&mut self, size: usize) -> Vec<Stereo<f32>> {

  // how much we need to prune
  let to_take = ((size as f64) * self.playback_rate) as usize;

  // to interleaved samples (original frame format is Vec<[f32, 2]>)
  let interleaved: Vec<f32> = self.take(to_take)
     .flat_map(|x| x.to_vec()).collect();

  // to complex DspVec data struct
  let mut complex = interleaved.to_complex_time_vec();

  // interpolation is HERE !
  complex.interpolatef(
    &mut self.interp_buffer,
    &self.sinc_function,
    (1.0 / self.playback_rate) as f32,
    0.0,
    16,
  );

  // re-frame the signal
  let chunked = complex.to_slice().chunks(2);
  let mut out = Vec::new();
  for chunk in chunked {
    let c = [chunk[0], chunk[1]];
    out.push(c);
  }

  // send full buffer
  return out;

}

Ah, that's interpolatef (it's in vector/src/vector_types/time_freq/interpolation.rs).

I'm not fluent in Rust, but I think it calculates the interpolation by:

Taking in a function (a real impulse response), and generating sampled vectors of that, with fractional offsets.
It uses that polyphase filterbank to do the arbitrary resampling

You could hence just "prepend" every chunk with conv_len samples, and then "throw away" the right amount of samples from the beginning of the result.

However, that might be hard: Since this function allows for arbitrary real resampling ratios, conv_length * interpolation_factor might not be an integer.

So, the easiest thing would probably to restrict yourself to rational resampling ratios, here, and always transport a number of samples $n \ge \text{conv_len}$ with $\mathbb N \ni n = \frac m{\text{conv_len}}$ so that you can always throw away an integer amount of samples.


Transform coding gain of two channel Haar Wavelet transform (first stage) Question

My interpretation

Work

I'm attempting to do this in MATLAB.

I know that for discrete Random Processes going into an LTI system the autocorrelation function of the output is $R_y[n] = h[n]\star h[-n]\star R_x[n]$ where $\star$ is convolution.
The two tap Haar filters are: h1 = [0.707 0.707] and  h2 = [0.707 -0.707] which are orthonormal to each other.

I also know that transform coding gain is defined this way:
$$G = \frac{1}{N}\frac{\displaystyle\sum_{i=0}^{N-1} \sigma_i^2}{\left(\displaystyle\prod_{i=0}^{N-1} \sigma_i^2\right)^{\frac 1N}}$$

where $N$ is the number of samples (3 in this case?) and $\sigma_i$ is the variance of the $i$th transform coefficient.

This is my code:

h1 = 1/sqrt(2).*[1 1];
h2 = 1/sqrt(2).*[1 -1];
y1 = conv(x,conv(h1,fliplr(h1)))
y1 = 4.6500   14.3000   19.3000   14.3000    4.6500
y2 = conv(x,conv(h2,fliplr(h2)))
y2 = -4.6500    4.3000    0.7000    4.3000   -4.6500

After decimation:

Y1 = 4.6500   19.3000    4.6500
Y2 = -4.6500    0.7000   -4.6500

Calculated coding gains:

GTC1 = 9.533/(417.3134)^(1/3) = 1.2757
GTC2 = -2.8667/(15.1358)^(1/3) = -1.1589

My questions:

Am I calculating $y_1$ and $y_2$ correctly? As I understand it these are supposed to be the autocorrelation functions for a WSS Gaussian random process (because WSS into LTI gives WSS out). After decimation, can i use $Y_1$ and $Y_2$ to calculate coding gain? I'm getting a negative coding gain for the second channel ($Y_2$), does that make sense?

EDIT: Just realized that the Haar filters need not be 2 tap filters. Now I am more lost than before.

I am not sure why you're trying to do this in MATLAB. This can be done quite easily with pen and papers.

Another way of seeing what you're doing here is the following:

Build vectors of size 2 from the sequence $X$
$$[X(0),~X(1)], [X(2),~X(3)], [X(4),~X(5)]\dots$$
Apply an orthogonal transform on each vector to obtain the corresponding transform coefficients $Y$.

The order 2 orthogonal transform you're using can mathematically be written as
$$\begin{bmatrix}Y_0 \\ Y_1\end{bmatrix} = \frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1 \\ 1 & -1\end{bmatrix}\begin{bmatrix}X(n) \\ X(n+1)\end{bmatrix}$$
where $Y_0$ and $Y_1$ are the transform coefficients and $X(n)$ and $X(n+1)$ are successive (hence correlated) samples from the process to be compressed. From this, you can easily compute $\sigma^2_0$ (the variance of $Y_0$) and $\sigma^2_1$ (the variance of $Y_1$). Indeed, because $X$ is zero-mean (and so are the transform coefficients), you know that
$$\sigma^2_i = R_{Y_i}(0) = E\{Y_iY^*_i\}$$
for $i = 1, 2$.
Then, you simply need to plug in the expression of $Y_i$ as a function of $X(n)$ and $X(n+1)$ to find $\sigma^2_i$ as a function of $R_X(0)$ and $R_X(\pm1)$. At the end, you should obtain
$$\sigma^2_0 = R_X(0) + R_X(1) = 19.3,$$
$$\sigma^2_1 = R_X(0) - R_X(1) = 0.7.$$
You can clearly see from this the energy compaction in the first transform coefficient, which will contributes to the coding gain. You now have everything you need to compute the coding gain (note that $N$ in your expression for the coding gain is the number of transform coefficients, 2 in this case).


digital IIR filter delay and oversampling Can oversampling decrease the delay of digital IIR filter? Imagine there is some digital signal going into processor that applies low pass filter.Lets say its 1 KHz sample rate and the filter is second order gaussian lowpass with -3db point at 100 Hz.

The putput of this digital filter will be delayed,this delay will vary with frequency.Can this delay,in the band of our interest which is 0 - 500 Hz for 1 KHz sample rate,be decreased if we increase the sample rate of the incoming signal?

What if the incoming signal was 2 KHz samplerate ( 2x oversampling ),would it decrease the digital filter delay?

I have read that the delay can be decreased by increasing bandwidth or by putting -3db cut off point higher in frequency.Is that true? I know making the low pass have higher frequency decrease the delay,but what about the bandwidth?

If we increase the sample rate,we increase the bandwidth,but in case of low pass filter it will be the bandwidth thats being attenuated.I have feeling that the author who wrote about the bandwidth thing meant the passband bandwidth and in that case its same thing as making the low pass higher.

I tried two filter simulating softwares,the Iowa Hills IIR and MicroModeler.com and got conflicting results.The Iowa Hills showed that the delay in the band of interest is increased slightly by oversampling while MicroModeler showed small decrease.

Does oversampling increase,decrease or not change the delay?

Note that a given digital filter impulse response $h[n]$ will have a corresponding equivalent analog frequency response $H(e^{j\omega})$ through the sampling relations; hence with the given sampling rate Fs.

This means that if you change the sampling rate of the input of this filter, then the effective analog filter response will also be changed, probably that's not what you want to have. So you have to redesign the digital filter to recreate the desired analog filter response.

In particular, when you increase the sampling rate by two, then the digital filter $h[n]$ should be redesigned such that its digital frequency response cutoff frequency is halved, so that the associated analog cutoff frequency would remain the same.

That means that, since the digital cutoff frequency is halved, the filter impulse response is doubled in length, the delay is doubled too. But don't worry, as the sampling period is also halved, the associated analog delay remains the same.

Of course, delay here refers to that associated with the decay time of the digital IIR impulse response $h[n]$.


Phase difference calculation for non divisible sample rate to frequency I am learning DSP and I am trying to simulate a simple phase difference calculator between two sinusoidal signals using GNU radio.

My flowchart looks like this:

As can be seen on the chart, one of the signals is using a delay block (controlled by the "delay sig" range GUI).

Other interesting points are:

fft_width variable:
Calculated by an expression - int(freq/20)
freq variable:
Controlled by a GUI range with steps of 1.2KHz. Min 1.2KHz, Max 9.6KHz
The skip head and Keep M in N blocks - used together to extract only the FFT bin in interest. To do so the skip head calculates the first N blocks to be skipped via this expression: int((samp_rate/fft_width))+2 . Idea taken from this video: https://youtu.be/GJKbD--MsLM?t=326

The problem:

This flowchart works when the samp_rate and freq are divisible - say samp_rate=10k and freq = 1k. However, for more of a real world scenario I have chosen variables that have lower ratio and are not divisible. In particular: samp_rate=11k and freq=1.2k. As a result of this the phase calculation fluctuates between 180 degree difference - for example 3.42 and 3.42-pi.

The question:

How do I make this flowchart work when the samp_rate and freq are not divisible? Furthermore how do I make it work for samp_rate of 11k, but different frequencies (chosen by the freq range GUI)?

What I tried:

Not much unfortunately, as I cannot find much information about this problem. That said, I read somewhere that clock synchronization via a PLL or Clock Recovery MM might help, but I do not know how to apply the relevant GNU radio blocks to this flowchart. Furthermore, if PLL is used would it not change the phase difference between the signals when trying to synchronize them?

Thank you for your time!

EDIT:
I tried putting a PLL ref output block between the throttle and FFT blocks:

For the parameters I used:

loop bandwidth: pi/200 (documentation states it should be either that or 2*pi/100)

min_freq & max_freq: (2*pi*1.1e3)/samp_rate (documentation states it needs to be in radians per sample)

My goal was to try to make it work for 1.2k frequency, but unfortunately I still face the same problems. How can I fix it?

Re-reference your FFT phase measurement reference point to the center of your data window by doing an FFTShift (N/2 vector rotate) before the FFT.  That moves the portion of data where there is a circular discontinuity (due to non-integer periodicity) away from the phase reference point of a typical FFT.


What is the difference between positive matrix coefficients and negative matrix coefficients of an audio? I have turned an audio file into an one dimensional array by using audioread function in matlab and found several positive and negative fractional coefficient values.What is the basic difference of this positive and negative sign in the matrix?

A sound wave is a pressure variation in a fluid (typically air) that varies with time and position. A microphone records this pressure variation: Positive means the local pressure is slightly higher than the static pressure and negative means it's slightly lower.


Signals and systems : why do we study causal signals? Till now I have read that causal signals are right sided and anti-causal, left sided.

Why did we need to classify a signal with respect to its position?
What is it's physical interpretation?
Similarly, why are LTI systems classified into causal and non-causal systems ?

An LTI system (or even a system that isn't L or TI) that is "causal" has a prayer of being realized in real time whereas an "acausal" system cannot ever be realized in real time because an acausal system responds to input from the future.  A causal system responds only to input of the present and/or input from the past.

For LTI systems, it can be shown that the output $y(t)$ depends strictly on the input $x(t)$ and the system impulse response $h(t)$.  This dependence is expressed as the operation we call "convolution":

$$y(t) = h(t) \circledast x(t) \triangleq \int\limits_{-\infty}^{+\infty} h(\tau) \, x(t-\tau) \, d\tau $$

or

$$y(t) = \int\limits_{-\infty}^{+\infty} x(\tau) \, h(t-\tau) \, d\tau $$

If and only if $h(t) = 0$ for all $t<0$, then the output $y(t)$ has no dependence on future values of the input, $x(t)$.

$$y(t) = \int\limits_{0}^{+\infty} h(\tau) \, x(t-\tau) \, d\tau $$

or

$$y(t) = \int\limits_{-\infty}^{t} x(\tau) \, h(t-\tau) \, d\tau $$

That is the salient difference between causal and acausal systems.

Now since the impulse response $h(t)$ is called "causal" if it is zero for all time before $t=0$ (meaning the impulse response does not respond before the impulse hits the system input), that definition was extended to all signals, even signals that are not an impulse response.

The reason that it is significant is that it turns out that the real and imaginary part of the frequency response of a causal system must be Hilbert transform pairs; they must be related to each other as Hilbert transforms of each other.


How to filter a regular ticking out of a digital (on/off) signal? I have a digital on/off signal sampled at a regular interval.  There was a problem with the sensor that produced the signal which has sometimes caused a regular on/off ticking when the sensor should be in the off state.  I.e. sometimes when the readings should be zero, there is a square wave instead (but with an unequal mark:space).

How would I go about designing a filter to remove this square wave (i.e. output a zero if there is this regular ticking occurring, output the input value otherwise)?

The ticking has a frequency of around 16 minutes, the events that the sensor was observing might last anything from one minute to a few hours.  I expect that this means that perfect recovery of the correct data is impossible, but I’d like to get as close as possible.

I would think some kind of notch filter is what I need, but everything I’ve found so far seems to be about filtering analogue signals.

(My memory of signal processing courses is a bit patchy, hopefully I’ve used the right terms here).

EDIT:

Thanks for your comments so far.  I've attached an image illustrating what I am trying to achieve.  Also I wrote the wrong frequency for the regular ticking - it's actually 16 minutes, not 8, but I doubt that affects the overall approach.

@Jan - yes, ideally I would just have a better sensor, but this is the data I have for now, I don't currently have the option of switching to a different sensor.

EDIT 2:

The data processing is happening offline, so non causal solutions are fine.

Are you processing this data offline or in real time? If you're doing offline processing, or can live with a significant delay, then you could apply a simple filter as follows:

From each input sample, check to see if the sensor reported "on".
If it did, search for a "few" 16-minute periods in each direction and see whether the sensor reported "on" there as well. This is a noncausal operation; if you're processing data in real time, you can search in both directions by first applying a large enough delay to the signal.
If there are "enough" consecutive "on" samples at the spurious period in either direction, suppress the "on" indication from the sensor.

This should suppress the unwanted square wave of "on" indications, but as indicated in the comments, it also could possibly suppress valid "on" reports. The "few" and "enough" parameters you'll have to select yourself based upon what is reasonable given the overall behavior of your system. Basically, you want to set "few" so that you're looking at enough spurious periods to make a good decision. You'll want to set "enough" small enough to suppress short bursts of the spurious period, but large enough to avoid suppressing too many valid "on" indications that happen to be separated by 16 minutes.


Sparsity Representation of a Signal Using the DCT Matrix I have a signal $\mathbf x$, and I need to know how to obtain the matrix which is the corresponding sparsity basis $\mathbf\Psi$ such that $\mathbf x = \mathbf{\Psi\theta}$, where $\mathbf\theta$ is sparse. I know that the suitable basis for my signal is the DCT basis. Given the signal, how can I obtain its DCT basis using MATLAB?

If you know that your ($N\times 1$) signal $\mathbf x$ is supposed to be sparse when DCT-transformed, then your $\mathbf\Psi$ should be the IDCT (inverse DCT). Then your $\mathbf\theta$ will be in the DCT ("frequency") domain and your $\mathbf x$ will be in the "time" domain.

You basically have two options for implementing this:

If you are using this in a context where it is OK to use the DCT as a function instead of an actual matrix, this may be more efficient, as your programming environment's implementation of the DCT is most likely more efficient in terms of storage and computational complexity than explicitly doing the matrix-vector product. In this case, look for the idct function;

x = idct(theta);

If you must implement it with an actual matrix-vector product, you can obtain the matrix $\mathbf\Psi$ by transforming an identity matrix of the appropriate dimensions:

Psi = idct(eye(N));
x = Psi * theta;

Expect that operations involving Psi * theta will be slower and more memory-consuming than idct(theta).


Entropy when having symbols with same probability I am trying to work out this problem:

Assume that we are transmitting data from a dictionary of 15 symbols. Let the
probability of sending one symbol be the same as that of sending any other
symbol. Calculate how often we need to fetch a symbol to transmit information at
a rate of 5 bits/s.

How can this be found?

HINT: Compute the entropy $H$ of your source and note that in your case it represents the  number of bits required to represent each symbol. Then from
$\quad H$ [bits/symbol] $\times$ $x$ [symbols/second] = $5$ [bits/second]
compute $x$.


Derivative of $l_1$ norm I want to compute the following derivative with respect to $n\times1$ vector $\mathbf x$.
$$g = \left\lVert \mathbf x - A \mathbf x  \right\rVert_1 $$

My work:

$$g = \left\lVert \mathbf x - A \mathbf x  \right\rVert_1  = \sum_{i=1}^{n} \lvert x_i - (A\mathbf x)_i\rvert = \sum_{i=1}^{n} \lvert x_i - A_i \cdot \mathbf x \rvert = \sum_{i=1}^{n} \lvert x_i - \sum_{j=1}^n a_{ij} x_j\rvert$$
So the $k$th element of derivative is:

$$\frac{\partial g}{\partial x_k} = \frac{\partial }{\partial x_k}\sum_{i=1}^n \lvert x_i - \sum_{j=1}^n a_{ij} x_j\rvert $$
$$= \frac{\partial }{\partial x_k}\bigg(\lvert x_1 - \sum_{j=1}^n a_{1j} x_j\rvert +\cdots+ \lvert x_k - \sum_{j=1}^n a_{kj} x_j\rvert + \cdots\lvert x_n - \sum_{j=1}^n a_{nj} x_j\rvert  \bigg)$$
$$ =-a_{1k}sign(x_1 - \sum_{j=1}^n a_{1j} x_j)-\cdots+(1-a_{kk})sign(x_k - \sum_{j=1}^n a_{kj} x_j)-\cdots -a_{nk}sign(x_n - \sum_{j=1}^n a_{nj} x_j)$$

And my questions:

Is this derivation correct?
How I can represent the answer compactly?
Can you introduce me a source to master this material?

Apart from a sign error, your result looks correct. The term with $(1-a_{1k})$ should have a positive sign. Also note that $\text{sgn}(x)$ as the derivative of $|x|$ is of course only valid for $x\neq 0$. If you take this into account, you can write the derivative in vector/matrix notation if you define $\text{sgn}(\mathbf{a})$ to be a vector with elements $\text{sgn}(a_i)$:

$$\nabla g=(\mathbf{I}-\mathbf{A}^T)\text{sgn}(\mathbf{x}-\mathbf{Ax})$$

where $\mathbf{I}$ is the $n\times n$ identity matrix.


Why make the phase of an FIR filter linear? I am new to DSP and filter design. It is easy to design an FIR filter with linear phase by making the coefficient sequence symmetric. However, why make the phase linear?

Because linear phase of a filter results in constant group delay according to the relation between group delay $\tau$ and phase response $\phi(\omega)$ of the filter:
$$ \tau = - \frac{ d\phi(\omega) }{ d\omega} $$

But then you will ask: why then constant group delay? And the short answer will be that constant group delay will preserve the relative positions (in time or space) of different frequency (narrow band) components of the signal to each other.

Many applications, such as image processing, may require such a fidelity on the processed signals. But not all, and to some degree, non-constant group delay such as provided by IIR filters are well tolerated in applications like audio...


what is a pneumatic signal exactly? I'm starting to learn about automatic process control, and there are types of signals transmitted in a control system. Among that is pneumatic signal that is a pressure signal. But I still don't know what a pressure signal means? How does it work? If you could give me a simple reference to learn more about signals(not a book to read, I'm already reading a book about control..)

For a pneumatic system, information is carried by the pressure of a gas (or air) in a pipe, instead of electrical current. In antique Greece,  pneuma was used in medicine to describe  air circulating for the  functioning of living systems.

I am not sure the term is very common, although there are several related concepts such as pneumatic actuator, pneumatic controller, pneumatic signal transmitter or amplifier, (fundamentals of) air logic, etc.

By a form of metonymy, I'd consider a pressure signal as the measure, as a given location  over time, in certain units (psi, psig for pound/square inch [gauge]).

Apparently, due to pressure/temperature conservation laws in constant volume, pneumatic controllers can also work on temperature.

In addition to the above link, the short document Pneumatic Systems describes in Chapter Different kinds of basic circuits: Flow amplification, Signal inversion,  Memory Function,  Delay function.


Time-variance and causality of $y(t)=\sum_{k=-\infty}^{+\infty}x(t)\delta(t-kT)$ Is the following system time invariant and/or causal?

  $$y(t)=\sum_{k=-\infty}^{+\infty}x(t)\delta(t-kT)$$

I think that it isn't time invariant and it is causal, and my thought is based on the fact that this is actually a sampling system that samples on $T-$intervals the input function. Is this correct and if yes, how can the time-variance be formally proven?

I agree with you on the system being a causal system since at a point in time $t$, the output is only made out of the input signal at the same point in time, weighted with an impulse (or not).

You can prove its time-variance by showing that a time-shifted input signal yields another output signal than the original input signal, for example:

$x(t) = \delta(t-\frac{T}{2})$ evaluates to $y(t)=0$ while
$x(t) = \delta(t)$ does not.


Identifying the correct source of packets I have a basic question regarding digital communication , especially pertaining to wireless internet .

A packet is basically a stream of bits . And a bit is detected as zero or 1 depending on whether the amplitude of a signal of certain frequency is crossing a threshold or not .
Now suppose a particular host or source $ A $ is sending this signals at a frequency $\ f $  to a receiver $ R $
Then it is possible that some device , $ B$ near the receiver can be placed to detect the frequency and regenerate a signal in that frequency and direct it towards the receiver .

In internet communication , I guess the packets will definitely be containing information about the source IP address, the IP of $ A $  but $ B $ can alter the bits so as to insert send its own IP  address in the packets .

In such situations , how do real receivers , identify the correct source ?

This is a very broad topic, and a full answer requires a few textbook chapters. A few pointers:

This problem is not solved in the physical layer; in other words, don't think in terms of frequencies, bits and packets.
The physical layer can handle some problems such as multipath, error correction and detection, and synchronization. The problem you're referring to is solved in the upper layers in the stack, and the solution depends on the particular network implementation.
For example, in WiFi there is an authentication protocol that ensures (with high probability) that every received packet comes from the actual source. A man-in-the-middle wouldn't be able to generate packets that would be recognized as valid by the receiver; again, this would happen somewhere in the network stack, not at the physical layer.


Real time signal processing use cases for eigenvalues of symmetric matrices I realize that this might be somewhat of an unusual and specific question.

I know that the eigenvalues of symmetric matrices are used in a number of ways in scientific computing, such as for finding numerical solutions to differential equations.

However, when considering real-time applications, I have only been able to find applications which also require the eigenvectors to be found. These include the MUSIC algorithm and principal component analysis.

Does anyone here know of any real-time applications/algorithms which only requires the eigenvalues of symmetric matrices?

EDIT:

The reason I am asking is that I have designed a hardware architecture for finding eigenvalues of symmetric matrices. I am looking for some kind of application where I can demonstrate its capabilities.

OK, you know that complex sinusoids are eigenfunctions (eigenvectors) of linear time-invariant (LTI) systems (which are linear operators). Guess what you see when you look at a spectrum: it's the eigenvalues corresponding to the eigenvector of any given frequency. So, basic signal theory uses eigenvectors to understand the fact that the frequency- and time-domain representation of signals are equivalent.

In MUSIC especially, the eigenvalues are of great importance – they signify the power of that eigenvector of the observed system (i.e. the power of the frequency component in case of that eigenvector spanning the signal subspace). So, MUSIC depends on the eigenvalues.

Basically, all symmetrical MIMO schemes that involve precoding decompose a matrix using the EVD into especially a diagonal matrix with the Eigenvalues as entries – and the question whether, and with what power, a decomposed "virtual" channel is used depends on the magnitude of the respective eigenvalue. So, MIMO, Precoding, and in the context of power assignment waterfilling for these channels depend on eigenvalues.

Now, for things that are "borderline" signal processing:

Graph theory has a relatively rich interaction with linear algebra. The Spectrum of a Graph is basically equivalently used to Spectra of Matrices.

Especially for random graphs, statements on bounds and distribution of eigenvalues of the (often unboundedly large) adjacency matrix are useful to describe degree distributions, and how reliably one can remove random edges from the graph without cutting the graph.

You can imagine the relevance of that to the theory of routing in communication networks; but in the case of random parity-check codes, the distribution of node degrees and cycle lengths is an important thing that can be at lest approximated through this.


DFT and integer valued basis functions'sf requencies In Matlab the function

W=dftmtx(N)

gives the DFT matrix of size N.

Each row is computed for an integer frequency k. $W_{k,n} = e^{-i2\pi kn/N}$, k-th frequency, nth sample if I am not mistaken.

I remember that I ve learned that the frequencies (k) are integers for the DFT. For the continuous Fourier transform it is a continuum of freq and times so everything they are real numbers both.

Could you develop on the usage of integer frequencies? and possibly relate it to the DCT where it seems to me that, e.g. the second basis vector is like half a cycle (as shown on my plot), and half a cycle is like a fractional frequency?

code for plot:

DCTm=dctmtx(8);
DFTm=dftmtx(8);
figure;plot(DCTm(2,:));hold on;plot(real(DFTm(2,:)))

Could you develop on the usage of integer frequencies?

When dealing with sampled discrete signals, there can only be an integer amount of resolvable frequencies within a captured interval.

When you take an $x(t), t \in \mathbb{R}$ and apply sample-and-hold on it at some sampling frequency $Fs$, it is turned into an $x \left[n \cdot \frac{1}{Fs} \right ], n \in \mathbb{N}$.

Theoretically, when $t \in \mathbb{R}$, you can "interrogate" any time instance or interval of $x$. But after sampling, when $n \in \mathbb{N}$, the smallest interval you can interrogate from your signal is $\frac{1}{Fs}$ or between successive samples $n, n+1$.

Even if you tried to sample faster after the sample-and-hold, the value you would get would be the last known value from the last time the sampling took place.

Therefore, within a block of $N$ samples, you get $N$ resolvable frequencies. Even if you tried to evaluate the Discrete Fourier Transform in more than the available $N$ frequencies, all that you would get would be interpolated values between the $N$ resolvable bins of it.

Which brings us to what Stanley Pawlukiewicz remarks: The fact that $n$ is integer, does not mean that you can only represent integer frequencies. Only that you can represent a fixed amount of them.

Just as you can only "see" $n$ and then $n+1$ discrete time instances of $x$, so you can only resolve some $f, f+ \frac{Fs}{N}$ discrete frequencies within it, where $N$ is the amount of samples you have collected.

  ...possibly relate it to the DCT where it seems to me that, e.g. the second basis vector is like half a cycle (as shown on my plot), and half a cycle is like a fractional frequency?

I hope that it is clear from the above discussion that "...half a cycle..." can well be N samples of some $f$ sampled at the right $Fs$ (?).

Hope this helps.


Symbol timing recovery in Python My objective is to demodulate QPSK signal. At the receiver I apply RRC filter and interpolate the signal to get the values of the signal (approximately) at the sampling instances. Then I am concerned with the symbol timing recovery to get the correct sampling timing. The filtered and interpolated signal is represented in the figure below.

Red vertical lines indicate the symbol intervals. Obviously, the symbol timing recovery is needed.
My idea is to apply PLL in order to track the phase error to zero and then sample the corrected signal at correct time instances.
So I have the following python code:
import pdb

class SimPLL(object):
    def __init__(self, lf_bandwidth):
        self.phase_out = 0.0
        self.freq_out = 20*pow(10, 6)  # Arbitrarily set, because I expect to get the input signal in the range of 20 MHz. Thus with this frequency I think the PLL can be locked faster
        self.vco = np.exp(1j*self.phase_out)
        self.phase_difference = 0.0
        self.bw = lf_bandwidth
        self.beta = np.sqrt(lf_bandwidth)

    def update_phase_estimate(self):
        self.vco = np.exp(1j*self.phase_out)

    def update_phase_difference(self, in_sig):
        self.phase_difference = np.angle(in_sig*np.conj(self.vco))

    def step(self, in_sig):
        # Takes an instantaneous sample of a signal and updates the PLL's inner state
        self.update_phase_difference(in_sig)
        self.freq_out += self.bw * self.phase_difference
        self.phase_out += self.beta * self.phase_difference + self.freq_out
        self.update_phase_estimate()

pll = SimPLL(f_symb//20) # I have read in one post, that with the bandwidth of the loop filter between f_symb/100...f_symb/20 PLL performs well
num_samples = 100
phi = 0
ref = []
out = []
diff = []
for i in range(0, num_samples - 1):
    in_sig = sig_interpolated[i]
    phi = np.angle(sig_interpolated[i])
    pll.step(sig_interpolated[i])
    ref.append(sig_interpolated[i])
    out.append(pll.vco)
    diff.append(pll.phase_difference)
#plt.plot(ref)
#plt.plot(ref)
#plt.plot(out)
#plt.plot(diff)

ref, = plt.plot(ref, label='ref')
out, = plt.plot(out, label='out')
diff, = plt.plot(diff, label='diff')

plt.legend(handles = [ref, out, diff], loc = 'upper right')

plt.show()

And the result I get is represented in the figure below:

So it seems like the PLL cannot lock the loop and the signal cannot be corrected.
Does anybody know, how to fix my code or can point me to the relevant example?
Thanks!

Phase offset and timing offset are two different things with two independent tracking loops to resolve. A Phase Lock loop eliminates carrier phase offset (and a frequency offset is a carrier offset changing with time). We can have phase offset with no timing offset, and we can have timing offset with no phase offset.
Once carrier offset is removed, I recommend creating an upsampled eye diagram to evaluate the signal quality and see if further timing recovery is all that is needed. The diagrams below specific to 16 QAM illustrate the difference between carrier offset and timing offset:
Here is the eye diagram of the real portion of the 16 QAM signal, where we see the four levels and we see the signal is sampled at 4 samples per symbol. The signal quality is excellent as given by the wide eye openings and small variation at the ideal sampling locations. However given the present sampling locations there is a clear timing offset where further timing recovery is needed. There is no phase offset in this example.

For comparison, after timing recovery the eye diagram looks like the following:

If we plot only the samples at the correct sampling locations (one sample per symbol) on an an IQ diagram to get the constellation, the above case including both real and imaginary values would be at the orange dots in the following graphic:

If there was a phase error in contrast to timing offset, the entire constellation would be rotated by the phase error, demonstrating phase error with no time error:

Note that if there is phase error even with no timing error the eye diagram will appear very corrupted still (often there is a phase rotation due to frequency offset which makes the eye unrecognizable, these must be corrected first before evaluating signal fidelity with an  eye diagram). Below shows the same eye as was done above for 16 QAM with no timing error, but with the phase offset as shown in the previous constellation:

A Gardner Loop is my recommended go-to for timing offset correction in that it can converge and remove timing errors even in the presence of very large carrier offsets. Then the carrier offset can be determined by simply measuring phase rotation from symbol to symbol using the correct sampling points as determined by the Gardner timing recovery loop.
More details on the Gardner Loop are provided here:
Gardner Timing Recovery for Repeated Symbols
And more details on Carrier Recovery are provided here:
High modulation index PSK - carrier recovery


A good textbook for designing signal filters Since a couple of months ago I started being actively involved in the area of control of dynamical systems.

In most cases, designing a controller for a given dynamic system will require the employment of digital signal processing techniques, especially in the area of signal filters.

Since I do not have a background in control engineering, I wonder if any of you may provide me with suggestions of sound textbooks dealing with signal filters into some detail.

The textbooks should ideally encompass:

Introduction to the most common filters in the area of DSP;
Explanation of their main characteristics in both frequency and time domain;
In which scenarios are they generally used (role or function of the filter).

Although being a naive question, I hope you may advise a couple of textbooks.

There are a lot of books out there, but if you are interested in Control and Signal Processing, I strongly suggest you take a look a Stephen Boyd Lectures from standford:

http://www.youtube.com/watch?v=bf1264iFr-w

There's the first one, the entire course is really valuable and he is a great Teacher.

Appart from That here's a good list of my preferred books on Signal Processing, some are more introductory, and some are more advanced:

Basics:

Signals And Systems - Oppenheim & Willsky
Digital Signal Processing - Proakis & Manolakis
Discrete Time Signal Processing - Oppenheim & Schafer

A bit more Advanced and Related to Control:

Optimum Signal Processing - Sophocles J. Orfanidis
Optimal Filtering - Brian D. O. Anderson & John B. Moore
Optimal Control: Linear Quadratic Methods - Brian D. O. Anderson & John B. Moore

Some on Statistical Processing (Advanced):

Statistical Signal Procfessing - Louis L. Scharf
Linear Estimation - Thomas Kailath, Ali H. Sayed & Babak Hassibi (Watch out this one is really tough)
Spectral Analysis of Signals - Petre Stoica & Randolph Moses

And some about Adaptive processing (Useful in every area of signal processing and control):

Adaptive Filter Theory - Simon Haykin
Adaptive Filters - Ali H. Sayed
Satistical And adaptive Signal Processing - Manolakis, Ingle & Kogon

And of course, the father of all Control Engineering Textbooks:

Linear Systems - Thomas Kailath

From your question i cannot really guess your background, but if you have never studied any topic like LTI systems, Digital Systems or Control, i would say you have a couple of years ahead of you to really learn and dominate the topic. The books i listed are the most referenced books in university courses, and in general the ones that everyone studies. They are for the most part, very theoretical, a lot of math and stuff like that. But it's a must have for reference. If you are looking for more instant satisfaction i would say try to google some cookbooks about digital filters or control, with the basics of PID controllers, linear regulators, etc.

Hope it helps


why is the signal divided into epochs in EEG classification? I'm new in EEG signal classification.
Studying the literature on this topic, I wonder why the EEG signals are divided in epochs, so, instead of classifying the whole signal all at once, we usually classify "windows" of signal, splitting the original one. I searched in the literature for a possible explanation, but I find nothing which answers to my question.

Basically, to perform a spectral analysis. Take a look to this: http://support.ircam.fr/docs/AudioSculpt/3.0/co/Window%20Size.html
I do not know if you refer to some specific approach, but usually, the main EEG classification approaches use frequency analysis (also combined with temporal analysis) to extract the features that will be used to train the classifier or infer the prediction.
Cheers


Increasing Recall rate for SIFT Is there a better way to increase Recall Rate when using SIFT features?
I am thinking a way to replace the NN1/NN2 ratio to account for slightly distorted objects.
Moving towards clustering and using BOW(Bag Of Words) seems a way but I need to do one-on-one match of objects in images rather than training and learning. This refrains me from thinking towards BOW.

Anybody got any idea?

With slight modification you might want to use RootSift:

http://www.robots.ox.ac.uk/~vgg/publications/2012/Arandjelovic12/arandjelovic12.pdf

Also the other steps in the paper will guide to improve the recall rate.

Cheers,


Calculate average mean FFT Magnitude in bins I am analyzing sounds of daily activities recorded by a smartphone. For example walking, getting up from bed, falling, running etc.. (one at the time).
Let’s take one them. My goal is to

convert from time domain to frequency domain
divide it in bins of 10 Hz wide
calculate the average FT magnitude in each of these bins.

My questions are:

Which FFT Magnitude spectrum should I use between the following graphs?
How should proceed to divide my signal (99144 total samples with 44100 sample rate) in 10 Hz bins? Can I achieve this just by using Matlab's pwelch function?

-> Here the frequencies of the single-sided magnitude are 0 – 22000 Hz -> Which means I should have 2200 bins of 10 Hz
-> Should I just create a new array and put the average mean of all values between 0 – 10 Hz,  10 – 20 Hz and so on?

[signal,fs]=audioread(filename); % fs = 44100
N = length(signal); % N = 94144

%% Single-sided magnitude spectrum with frequency axis in Hertz
% Each bin frequency is separated by fs/N Hertz.
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N; % conversion Hz - Bins
N_2 = ceil(N/2);
figure;
subplot(4,1,1);
plot(fax_Hz(1:N_2), X_mags(1:N_2))
xlabel('Frequency (Hz)')
ylabel('Magnitude');
title('Single-sided Magnitude spectrum (Hertz)');
axis tight

%% Single-sided magnitiude spectrum in decibels and Hertz
X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
subplot(4,1,2);
plot(fax_Hz(1:N_2), 10*log10(X_mags(1:N_2)))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)');
title('Single-sided Magnitude spectrum (Hertz)');
axis tight

%% Single-Sided Amplitude Spectrum
NFFT = 2^nextpow2(N);
Y_NFFT = fft(signal,NFFT)/N;
f_NFFT = fs/2*linspace(0,1,NFFT/2+1);
subplot(4,1,3);
plot(f_NFFT,2*abs(Y_NFFT(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

%% Pwelch defaut values
subplot(4,1,4);
pwelch(signal,[],[],[],fs);
title(sprintf('Pwelch defaut values'));

Here the answer for those facing similar problems:

if your Fs = 44100 Hz (in my case) and you want fft resolution of 10 Hz then use fft resolution N = 44100/10 = 4410.

You can then take just 4410 samples of signal or split the stream into segments (+ overlap) each = 4410 samples and directly average the bins (fft output).

FFT covers frequencies from -Fs/2 to +Fs/2 so you don't need to worry about single sided spectrum.

You better average the amplitude values from each bin directly(not the dB), then you can convert the mean to power in dB:

fft_mean = mean(ffts1,ffts2,ffts3,...);

plot(20*log10(f,abs(fft_mean))); %in dB


The Kernel of the Guided Image Filter I am trying to understand how to calculate kernel matrix for guided image filter. Following is the formula for kernel calculation.

where k is window withing pixels i and j belong i assumed that window k is a 3x3 matrix and take following I matrix to calculate kernel. omega is total number of pixels in window.
I am confused about index i and j so if I take a matrix 3x 3 that is as
$$
     I=   \begin{matrix}
        1 & 2 & 3 \\
        4 & 5 & 6 \\
        7 & 8 & 9 \\
        \end{matrix}
$$
now mean is 5 and variance is 6.67 taking epsilon as 0.002 and omega is 9 I calculated value for kernel for the first pixel so taking i = pixel(1,1) and all rest will be j including i== j case. then my output came as (I used excel sheet to do this I could have used MATLAB)

$$
     W=   \begin{matrix}
        0.111 & 0.111 & 0.111\\
        0.111 & 0.111 & 0.111 \\
        0.111 & 0.111 & 0.111 \\
        \end{matrix}
$$

when I changed values in matrix I output was still same so I doubt i understood something really wrong. Is there anyone who can explain me how to calculate W in this formula.
Most confusing term for me is that sum of k:(i,j) over window k. If anyone can explain me I guess I will be able to correct my calculation.

Reference

The paper - Guided Image Filtering (Available on Scribd).
Please check section 3.3 for kernel formula.

The Guided Image Filter is a filter which assumes Piece Wise Linear model of the image.
Basically, it tries to estimate the best Linear Estimator per patch / pixel of the image.

Since it is Linear its coefficients are based on the first and second moments of the local patch.

You can see the code in my Q42415 GitHub Repository which also includes an implementation in Kernel Form in ApplyLocalLinearKernel().
Here are some results (Output of Q42415):


2-PAM or 4-PAM if I want to send data as fast as possible? I would like to ask if I should use 2-PAM modulation or 4-PAM modulationif I want to send data as fast as possible given that in each case I have the same frequency spectrum available.

4-PAM transfers 2 bits ($log_2(4) = 2$) per symbol, while 2-PAM transfers 1 bit per symbol ($log_2(2) = 1$), so 4-PAM definitely transfers data faster.  The catch is that you need higher SNR to get the same bit-error rate (BER) that you would get with 2-PAM.


All-pass filter usage in automotive room correction This may not be the right place to ask this question but I'll give it a shot anyways. Maybe someone here is an audiophile.

I'm using a DSP called Audison Bit One HD which uses FIR filters for crossover and EQ. They recently included all-pass filters which allows me to adjust the phase of the crossover point. I'm wondering if anyone here knows how this will benefit for room correction.

I know what an all-pass filter does but I'm not sure how it could be useful with FIR filters which is a linear phase filter. Why would I want to adjust the phase of the crossover point for?

The manufacturer is no help.

Allpass filters are very helpful in automotive mainly for getting the imaging dialed in. You often need to time align the contributions of different driver groups (left, center , right) at a specific listening position. Due to the complicated acoustics of the cabin. the delays are often dependent on frequency and so all passes can be very helpful in correcting this.
Allpasses in cross overs are typically not that useful, since the range of overlapping frequencies is pretty narrow and you can typically get by with a single delay


Sampling at non-uniform intervals I measure noise in different systems by looking at voltage fluctuations. I wanted to use a Keithley SMU for that purpose, by storing the data in the buffer and then reading them out. However, I found that the data points are not measured at equally spaced points. While the difference is small, it is still visible. For example, the time intervals for a given trial where: 3.0 ms, 3.01 ms, 2.99 ms, etc, instead of a fixed 3 ms time interval.

I am wondering whether I can still overlook this error or not. If not, what is the main reason? Also, is there anything else I can do?

You still can do a lot. There is no perfect sampler, jitter always exists. Moreover, traditional systems include quantization which prevents (theoretically) the satisfaction of Shannon-Nyquist-Kotelnikov-Raabe... conditions.
If your system includes low-pass filtering, or/and the phenomenon you look at is very weak above 166 Hz, the variations are likely to be small enough for your purposes. I would propose to look at an harmonic behavior with the following three approaches:

Suppose the periods are even, then perform Fourier analysis,
Resample in a simple fashion at even samples (linear interpolation), then perform Fourier analysis,
Use a Fourier tool adapted to irregular sampling, like the Lomb-Scargle periodogram.

If you detect no meaningful difference in the features you extract, then you are fine (for the time being). You even can, routinely, use these methods in parallel to detect a weird behavior when they differ.

If you see meaningful difference then you might want to invest into irregular sampling processing methods. Luckily, in theory, when you have irregular samples, and a sufficient precision in the timing, you can extract information beyond Nyquist, at the cost of more involved methods. There is a handful of methods for Fourier or time-frequency methods, FIR results, etc.


Design of a Complex FIR Filter For apply least-squares linear-phase FIR filter design,with frequency domain specification is not symmetrical.

The pass-band error function,

$$E(\mathbf{h})_p=\int_{\omega_{p_1}}^{\omega_{p_2}}| \mathbf{c}^T(\omega)\cdot \mathbf{h}-D(\omega)|^2d\omega \tag{1}$$

The stop-band error function,

$$E(\mathbf{h})_s=\int_{\omega_{s_1}}^{\omega_{s_2}}|\mathbf{c}^T(\omega)\cdot \mathbf{h}|^2d\omega\tag{2}$$

h: unknown complex coefficients.

Total error function,

$$E(\mathbf{h})_t=w_1E(\mathbf{h})_p+w_2E(\mathbf{h})_s\tag{3}$$

$w_1,w_2$: weighting constants .

The equation $(1)$ : which represents a filter approximation by determine the error between actual and desired response in pass-band , why the method is different in stop-band (equation $(2)$)  ?

-In matlab ,does it the freqz function achieve the transpose operation ?

In the stopband(s), Eqs $(1)$ and $(2)$ are equivalent, because in the stopband the desired response equals zero: $D(\omega)=0$. The reason why you might want to split the error in passband and stopband error is to apply different weights, as shown in your Eq. $(3)$.

The function freqz computes the frequency response of a discrete-time filter on a grid of frequencies, given the filter coefficients. I'm not sure what you mean by "achieve the transpose operation".


Relation between $X(f)$ and $S_x(f)$ We know that for a signal $x(t)$, it is related to $R_x(\tau)$ as,  $$R_x(\tau) = \int_{-\infty}^{\infty}x(t)x(t-\tau)dt$$.
$\\$We also know that $$R_x(\tau) \rightleftharpoons S_x(f)$$
$\\$How do we mathematically relate $X(f)$ with $S_x(f)$, where $X(f) \leftrightharpoons x(t) $?

$$ \mathscr{F} \Big\{ x(t) \Big\} \triangleq X(f) \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2 \pi f t} \ \mathrm{d}t $$

$$ \mathscr{F}^{-1} \Big\{ X(f) \Big\} \triangleq x(t) = \int\limits_{-\infty}^{\infty} X(f) \, e^{j2 \pi f t} \ \mathrm{d}f $$

Assuming $x(t)$ is real (which means that $R_x(\tau)$ is also real).

$$\begin{align}
 S_x(f) &\triangleq \mathscr{F} \Big\{ R_x(\tau) \Big\} \\
\\
 &= \int\limits_{-\infty}^{\infty} R_x(\tau) \, e^{-j2 \pi f \tau} \ \mathrm{d}\tau \\
\\
 &= \int\limits_{-\infty}^{\infty} \left( \int_\limits{-\infty}^{\infty}x(t)x(t-\tau) \mathrm{d}t \right) \ e^{-j2 \pi f \tau} \ \mathrm{d}\tau \\
\\
 &= \int\limits_{-\infty}^{\infty} \int_\limits{-\infty}^{\infty}x(t)x(t-\tau) \ e^{-j2 \pi f \tau} \,\mathrm{d}t \,\mathrm{d}\tau \\
\\
 &= \int\limits_{-\infty}^{\infty} \int_\limits{-\infty}^{\infty}x(t)x(t-\tau) \ e^{-j2 \pi f \tau} \,\mathrm{d}\tau \,\mathrm{d}t \\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) \int_\limits{-\infty}^{\infty}x(t-\tau) \ e^{-j2 \pi f \tau} \,\mathrm{d}\tau \,\mathrm{d}t \\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) e^{-j2 \pi f t} \int_\limits{-\infty}^{\infty}x(t-\tau) \ e^{-j2 \pi f (\tau-t)} \,\mathrm{d}\tau \,\mathrm{d}t \qquad \text{let } u \triangleq \tau-t\\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) e^{-j2 \pi f t} \int_\limits{-\infty}^{\infty}x(-u) \ e^{-j2 \pi f u} \,\mathrm{d}u \,\mathrm{d}t \\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) e^{-j2 \pi f t}  \,\mathrm{d}t \int_\limits{-\infty}^{\infty}x(-u) \ e^{-j2 \pi f u} \,\mathrm{d}u \\
\\
 &= \mathscr{F} \Big\{ x(t) \Big\} \cdot \mathscr{F} \Big\{ x(-t) \Big\} \\
\\
 &= X(f) \cdot X(-f) \\
\\
 &= X(f) \cdot \big( X(f) \big)^*  \\
\\
 &= \Big|X(f)\Big|^2  \\
\end{align}$$


How to get a non-equally spaced FFT back into the time domain I have a signal that I STFT and then filter using an ERB spaced filterbank. At some point after this I want to get the signal back into the time domain, how can I go about this? Using a standard iSTFT function won't work because it assumed linearly spaced frequency bins, AFAIK? I've put a code snippet below.
I'm also not sure what what to tag this question as apart from fourier-transform
Y = stft(sig) # Y.shape = (1025,4000)
fb = filterbank() # fb.shape = (20,1025)

Y_erb = matrix_multiply(fb,Y) # Y_erb.shape = (20,4000)
```

The short answer is that you can solve a least squares problem with the input signal as the decision variables.
About Nonequispaced FFT
There is a related question in stack overflow.
If you are not using python, you can take a look into the nfft github.
Here will find some information. Section 4.6 covers the inversion of such transforms.


Seed Region Growing Threshold I have this matrix

The following points

R1 = (1,7) = 0
R2 = (3,2) = 6
R3 = (5,5) = 3
R4 = (7,3) = 1

And using this equation |g(seedi) – g(pixel)| < T
for T=1

And this is the solution for threshold = 1

I don't understand why all the ones for example, belong to R1 = α

Based on the equation we have |0 - 1| < 1, which means it shouldn't become α.

Given your seed locations and the threshold of 1, that means that not every pixel can be labelled. See below.
$$
\begin{array}{ccccccc}
5 & 6 & 1 & 1 & 1 & 0 & \require{enclose}\enclose{circle}[mathcolor="red"]{\color{black}{0}}\\
6 & 7 & 7 & 1 & 1 & 1 & 6\\
3 & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & 7 & 7 & 7 & 6 & 5\\
4 & 3 & 7 & 7& 7 & 5 & 3\\
3 & 3 & 2 & 1 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 3 & 4\\
2 & 2 & 1 & 2 & 1 & 2 & 3\\
1 & 1 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 1 & 1 & 2 & 2\\
\end{array}\\
\Huge \updownarrow\\
\begin{array}{ccccccc}
5 & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} &\enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \require{enclose}\enclose{circle}[mathcolor="red"]{\color{black}{0}} & \require{enclose}\enclose{circle}[mathcolor="red"]{\color{black}{0}}\\
\enclose{circle}[mathcolor="green"]{\color{black}{6}} & 7 & 7 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="green"]{\color{black}{6}}\\
\enclose{circle}[mathcolor="blue"]{\color{black}{3}} & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & 7 & 7 & 7 & \enclose{circle}[mathcolor="green"]{\color{black}{6}} & 5\\
4 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 7 & 7& 7 & 5 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}}\\
\enclose{circle}[mathcolor="blue"]{\color{black}{3}} & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 2 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & \enclose{circle}[mathcolor="blue"]{\color{black}{3}} & 4\\
2 & 2 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 2 & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 2 & \enclose{circle}[mathcolor="blue"]{\color{black}{3}}\\
\enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & \enclose{circle}[mathcolor="cyan"]{\color{black}{1}} & 2 & 2\\
\end{array}\\
$$
It looks like there should be 4 labels, not three... and, yes, it looks like the 0 and 1 seeded regions have been merged.

In addition to the lack of labels for all pixels (3 vs 4), I wonder if it's because the calculations are being done in (double precision) floating point so that $|1.00 - 0.00| < 1.00$