What are the problems with designing an FIR filter using FFT? I'm trying to get an understanding of the relationship between an FIR filter designed from "first principles" using a filter kernel with convolution, and a filter designed in one of two ways using FFT (see below).

As far as I understand, the impulse response of an FIR filter is the same thing as the filter's convolution kernel. (Correct me if I'm wrong.)

Also, in my understanding, the component frequencies (i.e. Fourier transform) of an FIR filter's impulse response is the same thing as the filter's frequency response. And, therefore, the inverse fourier transform will give me back the impulse response.(Again, correct me if I'm wrong).

This leads me to two conclusions (ignoring phase response, or assuming linear phase response):

I should be able to design an FIR filter of arbitrary frequency response by "drawing" my desired frequency response, taking an IFFT to get the impulse response, and using that as my convolution kernel.
Alternatively, I should be able to create a filter by taking the FFT of the input signal, multiplying by my desired arbitrary frequency response in the frequency domain, and taking an IFFT of the result to produce the output signal.

Intuitively, it feels like 1 & 2 are equivalent, but I'm not sure if I could prove that.

It seems like people (and DSP literature) go to great lengths to design FIR kernels with predefined responses, using complicated (to me) algorithms like Chebyshev or Remez (I'm throwing out some names I've read, without really understanding them).

Why go to these lengths, when an FFT/IFFT transform exists for every possible FIR kernel?
Why not simply draw the exact frequency response you desire, take an IFFT, and there's your FIR kernel (method 1 above)?

One reason you see people designing FIR filters, rather than taking a direct approach (like both 1 and 2) is that the direct approach usually fails to take into account the periodicity in the frequency domain, and the fact that convolution implemented using an FFT is circular convolution.

What does this mean?

Suppose you have a signal $x = [ 1, 2, 3, 4]$ and a filter impulse response (convolution kernel; you are correct they are the same) $h = [ 1, 1 ]$.

The convolution $y = x *h$ is $[1, 3, 5, 7, 4]$, a 5-length vector.  If you use the FFT (of the wrong length, 4) then the answer you get is $[3, 5, 7, 5]$.  The reason for the difference is that the result of linear convolution of these two is length 5, but the result of circular convolution is whatever the FFT length was.

If the FFT length is greater than or equal to the length of the result of linear convolution, then the two are the same.  Otherwise, the two are not the same (unless the data somehow conspires to make it so, e.g. if one signal was zero).


I/Q Sampling and Demodulation A bit stream is mapped to symbols during baseband modulation, each symbol (say for QPSK) represents two bits. Each symbol has a corresponding I & Q value, which in QPSK, are used as the amplitude for a cosine and sine carrier wave, which are then summed for the QPSK passband signal. Where does the baud rate come into play? If it is the rate at which the bit stream is processed and modulating the carrier, how do I account for baud rate if I just have a bit vector converted into a complex vector of symbols, in say, Matlab?

When I have a receiver digitizing the passband signal I receive into I/Q samples, do these values correspond to the original I/Q values representing the symbol? If the sample rate 4x the baud rate, I should receive 4:1 received samples:transmitted symbols? I would then down sample into 1:1 Rx:Tx symbols, then map the symbols back to bits using the Tx constellation (assuming no rotation)?

The baud rate/symbol rate is the inverse of the time distance $T_\mathrm{s}$ of two consecutive symbols. It determines the bandwidth of the modulated signal and thus comes into play when simulating frequency dependent components like filters, a wireless channel, a digital-to-analog converter, and so on. If in a simulation system the carrier modulation is simulated in the passband the ratio of baud rate and carrier frequency has to be taken into account. Imagine your carrier is modelled by $\sin(2\pi (f_\mathrm{c}/f_\mathrm{s}) n)$ with a period of $N=f_\mathrm{s}/f_\mathrm{c}$, where $f_\mathrm{c}$ is the carrier frequency and $f_\mathrm{s}$ is the sample rate of the simulation system. In the simulation system, one symbol then has to have the length $N_\mathrm{s}=T_\mathrm{s}f_\mathrm{s}$ or $N_\mathrm{s}=R_\mathrm{s}f_\mathrm{s}$, where $R_\mathrm{s}$ is the baud rate.

If however, you "just have a bit vector converted into a complex vector of symbols" without further upsampling, then the baud rate does not have to be taken into account.

Usually, communication systems are simulated in the equivalent baseband, i.e. the generally relatively fast oscillating carrier is not modelled by a sine wave. Instead, the signal is thought of being shifted to the carrier frequency and only the frequency band of interest that contains the modulated signal is taken into account. Here, the baud rate again determines the width of this band.

If the receiver applies oversampling, i.e. if the sampling rate $f_\mathrm{A}>R_\mathrm{s}$, then the received signal has to be downsampled to reconstruct the transmitted symbols. Oversampling is usually applied for clock recovery or a fractionally spaced equalizer. In a simulation system you could assume perfect sample time and $f_\mathrm{A}=R_\mathrm{s}$. But this depends on what should be simulated, of course.

Update

Generally, in a receiver applying oversampling, the received samples corresponding to a symbol are not identical due to various impairments: noise, inter-symbol-interference, impulse shape, etc. To reconstruct the transmitted symbols, several steps are necessary including clock recovery, synchronization, channel estimation and equalization, and resampling.


Power Spectral Density from Probability Density Function The samples of a signal $x[n]$ are i.i.d. and follow a triangular pdf with $a = 0,\ b = 2,\ c = 1$:

The DC-power of the signal is
$$\mu_x^2 = \big(\mathbb{E}(X)\big)^2 = \left(\int_{-\infty}^{\infty} x f_x(x)dx\right)^2 = 1$$
the total power of the signal
$$\mathbb{E}(X^2) = \int_{-\infty}^{\infty} x^2 f_x(x)dx = \frac{7}{6}$$
and the variance (AC-power) as
$$\sigma^2 = \mathbb{E}\big((X-\mu_x)^2\big)= \mathbb{E}\left(X^2\right)-\mu_x^2 = \frac{1}{6}$$.
Thus I thought, that the PSD of the signal has the following form:
$$ S_{xx}(e^{j\omega}) = 2\pi \delta(\omega) + \frac{1}{6}$$
as the AC-power results in a constant value over the whole spectrum and the DC power in a Dirac-Delta at $\omega = 0$. Furthermore, the integration over the PSD should return the total power of the signal:
$$ \frac{1}{2\pi} \int_{-\pi}^{\pi} S_{xx} dx = \mathbb{E}(X^2) $$
Is my solution correct or did I miss something? Can this procedure be done for every probability density function under the assumption that all samples of $x[n]$ are iid?

If the $x[n]$ are i.i.d. random variables with mean $\mu$ and variance $\sigma^2$, then the power spectral density of this discrete-time random process (which is effectively white noise plus a possibly nonzero mean $\mu$) is what you have calculated.  The shape of the common pdf of the random variables is irrelevant in all this except insofar as the shape determines the mean and variance; for example. you would get the same PSD if the $x[n]$ were Gaussian random variables with mean $\mu$ and variance $\sigma^2$.


Estimate Variance of Additive White Gaussian Noise (AWGN) Given Multiple Realizations with Different Mean Given $ N $ images of the same scene, where each image is corrupted by additive white Gaussian noise of the same variance $ {\sigma}^{2} $. How can $ {\sigma}^{2} $ be estimated?

So we have an Image $ X \in \mathbb{R}^{m \times n} $ which we have few realizations noisy variations of it  $ {\left\{ {X}_{i} \right\}}_{i = 1}^{N} $:

$$ {X}_{j} = X + {V}_{j} $$

Where $ V \in \mathbb{R}^{m \times n} $ is AWGN such that $ {V}_{j} \left( i , j \right) \sim \mathcal{N} \left( 0, {\sigma}^{2} \right). $

Ex., for $N=2$, let $y$ be the ground truth scene and $x_n$ be the noisy relizations, then $x_1=y+ν_1$ and $x_2=y+ν_2$ where $ ν_1, ν_2 \sim N(0,σ^2) $ and hence $x_1−x_2=ν_1−ν_2:=ν$ where $ν∼N(0,2σ^2)$. Therefor $σ^2$ can be estimated by calculating the variance of the difference of the two realizations, and then dividing by two.

Basically each pixel is a realization so all you need is to work in the 3rd dimension (Though you can also get better by using the Spatial Data).

Method 1

So the trick here is to use the multiple images to estimate the Mean (True value) of each pixel and then calculate the STD on all samples (numRows * numCols * numRealizations).

Assuming we have single channel image. So we pack all give images into a tensor tI with dimensions: numRows * numCols * numRealizations.

We calculate the mean per pixel (Averaging on the 3rd dimension) and then subtract each image from the calculated average image.
Then we're left with many realizations of the noise (Well, noise and the left over from the estimated mean error).

Method 2

Estimate the STD in the 3rd dimension per pixel. Them average all numRows * numPixels estimations.

Method 3

Just as method 2. But since the property which obeys to linear operatiosn is the Variance calculate the variance along the 3rd dimension, average it over all pixels and then take the sqrt() of the average Variance.

Summary

Here is a graph showing the estimated noise STD as a function of the number of realizations:

As can be seen, for an image with many pixels (Say more than 20,000 or so) Method 3 becomes almost perfect.

The full code is available on my StackExchange Signal Processing Q61273 GitHub Repository.

Remarks

Method 1 - Basically the Maximum Likelihood Estimator for the problem. It is biased as this is the property of the ML for the Variance / STD.
Method 2 - Works nice, really intuitive.
Method 3 - Basically the MLE with the non biased estimation for the variance. As can be seen. It is also the best in real world.

Yep, Method 3 is much better only because it divides the data by N while Method 1 divides by N - 1. But it has significance impact on performance.


identification of sounds within audio signals I am looking for methods to identify sounds within audio files.  I was hoping to find some ready-made user-friendly software to do this but haven't found it yet.  I am familiar with Matlab and some C++ and have mech. eng. background so I can do some coding work myself to implement a solution but I'm not a serious programmer.

From my online searches I seem to have found there are 2 potential methods; cross correlation and fingerprinting.  It seems to me that the audio fingerprinting method is potentially more accurate than raw cross correlation.  At this point this investigation is for my own hobby purposes to be implemented with a robot.

Can anyone with related experience comment which approach may be more successful?    It seems there are several audio fingerprinting methods and some open source codes even for matlab so I am intrigued, however, these methods are mostly developed to identify entire songs.

The sounds I want to identify are short (from < 1s to ~5s).  I was thinking if the analysis sample rate was increased to fill the short sound file that maybe the method could work with very short files such as these sounds, but I do not know the implications for the rest of the techniques of fingerprinting.  Does this seem like a reasonable / possible approach to anyone?  If not are there references for other potentially applicable DSP techniques for this I can research?

Thanks.

Audio Fingerprinting can be used to match near identical recordings of sounds with a bit of additive noise. It will not work for similar sounds, cover versions etc.

Without knowing your application, I'd recommend giving Dan Ellis' audio fingerprinting code a try. You'll find MATLAB and Python versions here:

http://labrosa.ee.columbia.edu/matlab/audfprint/


Seeking compressive sensing imaging demo in MATLAB I am looking to experiment in MATLAB with image restoration/reconstruction using compressed sensing. I am relatively new to compressive sensing, and was inspired by the following Wired article to learn more about it:

Fill in the Blanks: Using Math to Turn Lo-Res Datasets Into Hi-Res Samples

In the Wired article linked above, the authors show a low-res picture of US President Obama, and snapshots of the image as it progresses through various stages of "recovery," ultimately ending up with a clean-looking image.

I have read quite a few papers, and even have coded a simple algorithm for sparse signal recovery in MATLAB, but I would like to be able to "play around" with a low-quality image in MATLAB, using some compressed sensing algorithms. I think it would help solidify my understanding of what is going on by incorporating a visual example.

So does anyone know of any such MATLAB demos/examples which provide the MATLAB code and a sample image (or images) that would permit one to experiment with the recovery process?

You can start with the Compressive Imaging Code code by J. Romberg, illustrating the paper  "Imaging via Compressive Sampling".

Another great source of information and codes is on Nuit Blanche.


Slow and fast filter According to my knowledge a filter with low cut off frequency is called slow and with hgih cut off frequency is called fast filter.
I am trying to understand more about these two. I know what happens to signal when it goes through these 2 filters(separately). But I wanted to understand more about these filters. I simulated a 1st order PT1 filter for both very low fc and very high fc.
Below are the bode diagram for slow and fast filter respectively.
slow:

fast:

I'd be greatful if any of you would give me what's actually going on and the difference between these two.

In the simplest terms, the slowness or fastness of a filter refers to its time domain characterization. Specifically a fast filter will have a short impulse response so that the output goes to steady state as quick as possible in time, whereas a slow filter will have a long impulse response, whose steady state output occurs after a longer delay, so that it will not be responsive to transient changes in the input.

Now according to the time-bandwidth product feature of LTI filters (which is based on a Fourier transform property) if a filter has wide frequency support, then its time support will be narrow, and vice versa. As a consequence, low cutoff frequency filters will have narrow frequency support and hence large time support. But a large time support is equivalent to long impulse response which is stated to be a slow filter as above.

Therefore a low cutoff filter (a narrow frequency support) will be a slow responding one, whereas a larger cutoff filter will be fast responding one.


FM demodulation amplifier with saturation I'm trying to understand how FM demodulation works. I read a lot on the internet and as far as I understand, the main process is to extract the phase difference (that is also the output) between a quadrature signal and the input signal, and that use this difference to keep those two signals locked in phase.

However the output signal amplitude depends a lot on the input signal strenght. As I understand, the input of the phase locked loop should be constant in amplitude.

Besides on this article on wikipedia it's written that FM recivers do not need a AGC, but it's enough to use a "high-gain IF amplifier which is intentionally driven into saturation". How can this be right? I mean.. If the amplifier is driven into saturation, the input wave will look more like a square wave and all the quadrature equations won't be true anymore! Can anybody explain me this step?

EDIT 1

I think that the input should be constant in amplitute because(I'm thinking about a PLL):

$ FM_in(t) = A sin( W_c * t + \phi(t) ) \\
 FM_in(t) * cos(W_c*t) = {A \over 2}(sen(\phi(t)) +  cos(2W_c + \phi(t)) ) $

So filtering out the high frequency cosine, I get:

$ DEM_Out(t) = {A \over 2} sen(\phi(t)) \simeq {A \over 2}\phi(t) $

And so far so good, but to keep in lock the loop I should feed this value to the VCO, and if A doubles(as I understand), the phase shift of the VCO doubles, too so that the loop won't be in lock anymore. I hope I have explained what I mean!

The receiver's output does not depend on the input amplitude.  From the point of view of the receiver, the output only depends on the frequency/phase of the input.

An fm signal in an ideal environment has a constant amplitude.  Another way of putting it is that fm modulation has a constant envelope.  Often this is advantageous because it allows the transmit amplifier to run at saturated output and not impact the modulation since amplitude is not varied as opposed to QAM and other non-constant envelope modulation types.  The same should hold true for the receiving side of the equation.  In practice I imagine you'd still need to watch out for overloading the receiver and causing harmonics and other nasties.

A write-up on my pyFmRadio blog post illustrates FM's amplitude independence from the modulating audio source.


Reconstructing the song using its dominant frequencies I have found out 20 dominant frequencies of a song by splitting it into 20 sound-clips and applying Fourier transform on each one of the sound-clips. Now I am trying to reconstruct the song back from these 20 dominant frequencies.

info_ = audioinfo('../wws.wav'); % 11 secs long song

Fs = info_.SampleRate;
chunkDuration = 0.6;
numSamplesPerChunk = chunkDuration*Fs;

freq_vec = [];
for startLoc = 1:numSamplesPerChunk:info_.TotalSamples
    endLoc = min(startLoc + numSamplesPerChunk - 1, info_.TotalSamples);
    [x, Fs] = audioread('../wws.wav', [startLoc endLoc]);
    x = x(:,1);
    dt = 1/Fs;
    t = 0:dt:(length(x)*dt)-dt; % time scale
    L = length(t);
    nfft = 2^nextpow2(L);
    y = fft(x, nfft);
    y = abs(y.^2);
    y = y(1:1+nfft/2);
    [v,k] = max(y);
    f_scale = (0:nfft/2)* Fs/nfft;
    freq_vec = [freq_vec, f_scale(k)];
end

freq_vec is the vector of dominant frequencies of the song. My question is how can I use these frequencies to build back the song.

Short answer: You can't.

Sound is more than its dominant frequencies, and you threw away anything but that.


Removing drift from integration of accelerometer data I am trying to get a positional data from the accelerometer data using the following steps:

Re-zero the accelerometer value
Removing mean from accelerometer value
First integration of accelerometer value to get velocity
Removing mean from velocity
Integrate velocity to get position

, and the following is the result from every integration value:

Can anyone suggest on how to remove the drift appear on the velocity value? Some suggest to push the velocity value to zero when accelerometer value is near to zero but I don't think this is very robust.

Double integration amplifies any offsets, non-linearities and noise.  These can't be removed without the use of some type of external reference point measurements (e.g. not from the accelerometer) or other information.  Calibration might be able to remove most of the ones that do not vary with time (temperature, etc.)

For instance, if the size of the room is known, you can robustly clamp position to be inside the walls and above the floor.


Speckle pattern - speckles detection I have the following problem. I have a speckle pattern (of light) like on the image below (left).
 I have to locate as precisely as possible the centers of all speckles.

That would not be a problem if those speckles would have some constant, known shape, e.g. Gaussians.

But they are irregular and I do not really know how to do this.

What I do now is:

I look for the highest point in the image and Im registering coordinates of this point.
I mask (NaN) the circular area near this point - more or less of the size of a typical speckle (image, 2nd subplot) and
I'm searching for another highest point. If there are no NaNs in the circular ares of this point, I register its position also. And so on, hundreds of times.

But it is not working to well. Is there any standard method in signal analysis for localization of individual speckles in the speckle pattern?

Image:

Left - speckle pattern
2nd  - masked areas
3rd  - registered positions if the masks can overlap
4rd  - ... if they cannot.

Essentially, isn't this just looking for local maxima in the image? Have a look at the following code. It dilates the image, i.e. putting each pixel the value of the maximum of its neighborhood. Then, it checks which pixels actaully were the maximum and marks them with a circle:

import skimage
import skimage.io
import skimage.morphology
import numpy as np
import skimage.draw
import scipy.ndimage

I = skimage.io.imread("Capture.PNG", as_grey=True)

selem = np.ones((9,9))
D = skimage.morphology.dilation(I, selem)
plt.figure(figsize=(12,12))

plt.subplot(221); plt.imshow(I); plt.title("Original image")
plt.subplot(222); plt.imshow(D); plt.title("dilated image")

maxima = I==D
plt.subplot(223); plt.imshow(maxima); plt.title("Detected Maxima regions")

labels, nlabels = scipy.ndimage.label(maxima)

plt.subplot(224); plt.imshow(labels)
print (nlabels)

C = scipy.ndimage.measurements.center_of_mass(maxima, labels, 1+np.arange(nlabels-1))
I2 = I.copy()
for x, y in C:
    rr, cc = skimage.draw.circle_perimeter(int(x), int(y), 5, shape=I2.shape)
    I2[rr,cc] = 1
plt.imshow(I2); plt.title("Center of mass of maxima"); plt.colorbar();

what is missing is e.g. a threshold to ignore small maxima.


The higher the sidelobes of a signal frequency spectrum, the greater the "cross-talk" between transmission channels Why is that 'The higher the sidelobes, the greater the "cross-talk" between channels' when transmitting signals via hardware electronics?

where,

Cross-talk:  a signal transmitted on one circuit or channel
of a transmission system creates an undesired effect and/or interference in another
circuit or channel.
Sidelobes: The sidelobes (spectral leakage) introduced to signal frequency spectrum after
application of window functions to the transmission signal.

Convolution theorem states convolution in time domain is the same as multiplication in Fourier domain, also the multiplication in time domain is a convolution in Fourier domain.

So when you truncate a signal (multiply it by a window), your signal's spectrum will be convoluted with the window's spectrum which widen your signal's spectrum depending on the main lobe width and also replicate signal spectrum depending on the side lobe height which is called frequency leak.
if you consider different channels over different bandwidth, when you truncate the signal the frequency leak will cause channel cross-talk.


Phase invariant detection of low frequency spikes I am a physics graduate student studying a new physical device via simulations. Thanks to SE DSP, reading Lyon's Understanding DSP as suggested (not all the way through it yet though), and a bunch of work, I have been learning, and enjoying, DSP. I’m at a point where I am completely stumped by something, and need some help.

For this project, signal processing appears very simple, as seen in this diagram:

However, what happens inside the "Physics Device" somewhat complicated. Simply put, when the external signal $f(t)$ and the physics device interact, there is some destructive interference. To me, this means there will be a narrow low frequency spike at the time where there is interference. When there is no destructive interference, there is a lot of stuff going on. The outputs relating to the physical device are very ugly compared to a communication system; there are many harmonics and other things happening.

When the physical device I'm simulating is running, and an external microwave signal $f(t)=A\cos(\omega_1 t +\phi_1)+B\cos(\omega_2 t + \phi_1)$ is introduced, destructive interference creates some low frequency oscillations at particular locations:

Input #1: $f(t)=A\cos(\omega_1 t +\phi_1)+B\cos(\omega_2 t + \phi_1)$

Output #1:

Th upper graph shows raw data, and the lower shows this after I apply a Low pass filter. If I use a low-pass FIR filter near DC, it gives me nice peaks to identify where the low frequency oscillations occur.

In these, the first peak near $60\textrm{ ns}$ is related to the first term, and the second smaller peak near $100\textrm{ ns}$ is related to the second term. Also, $B=0.2A$, so the second peak should be shorter than the first peak.

If instead of phase shift $\phi_1$ in the first term, the signal has a phase shift of $\phi_2$, or as: $f(t)=A\cos(\omega t + \phi_2) +B\cos(\omega_2 t + \phi_1)$, we get similar looking raw data:

Input #2: $f(t)=A\cos(\omega_1 t +\phi_2)+B\cos(\omega_2 t + \phi_1)$

Output #2:

But after the low pass filter, the peaks are not-so-beautiful.

Is there an algorithm that will allow me to determine when an external signal is introduced into the physical system regardless of the initial phase shift of the external signal? My eyes can identify from the raw data where the peaks should be, but a low pass filter doesn’t do the job consistently.

I have tried using signal energy, as $\int_t^{t+T} \lvert x(t’)\rvert^2 dt’$, but this fails in two regards.

It does not see the second peak
And it requires some threshold, which doesn’t work well as $A$ and $B$ can differ by several orders of magnitude.

Is there a good algorithm to get nice thin peaks from both signals? Even more, is it something that can be built with existing technology?

I can post data for the above graphs. How can I post it? They are small lists, with ~ 100,000 elements.

I would use a band-pass filter centred at the frequency of the peaks (10ns period?). That should give you more of a peak in the second signal.
Use a quadrature filter pair (band-pass filter and its Hilbert transform) to gain phase invariance. The output is called the analytic signal. Its amplitude is a local phase-invariant energy measure.

In MATLAB this would be

f_a = hilbert(conv(f,g,'same'));
amplitude = abs(f_a);
phase = angle(f_a);

where f is the signal and g is the band-pass filter.


Assumption of SINR measurement I am assuming there is an Access Point (AP) and IoT Nodes in NOMA uplink scenario. Following is assumed:

the AP does not know the channel gain of IoT Nodes
the transmission power of nodes is also unknown to AP

How can we measure the SINR of the nodes: as we do not have the Power x Gain product without the knowledge of Power and Gain of IoT nodes.
Is it right to assume that the receiver of the AP won't be able to get any legible data (information) unless the SINR is above a certain threshold?

How can we measure the SINR of the nodes: as we do not have the Power x Gain product without the knowledge of Power and Gain of IoT nodes.

Not at all, you have zero information about the N in SINR, and without channel knowledge you also can't have any information on the S, nor the I, even if you knew the transmit powers (which you don't).

  Is it right to assume that the receiver of the AP won't be able to get any legible data (information) unless the SINR is above a certain threshold?

Yes. If you're into NOMA, you should be aware of a bit of basic information theory (otherwise, it's hard to understand why you'd even ever want to go non-orthogonal):

In this case, a very basic corner case of the Shannon Limit holds true. Some call it the ultimate Shannon Limit, but what it really just says is that if your bit energy is below -1.59 dB of the energy of your redundancy, you can't get any information from your reception.

In a multi-user system, one can (often iteratively) gain information about the interference¹, but if noise² alone pushes the SINR below -1.59 dB, you've lost and better redesign your system.

 ¹Thus reducing restoring the receiver's ability to "remove" parts of the entropy the interferer adds to the received signal.

²"Noise" actually being all entropy in the reception that doesn't come from the information source. That incorporates all entropy of interference that you can't "cancel", however that happens.


Plotted ECG signals are not around Amplitude 0 line For my project I have to have few datasets of ECG signals and I did get them from Physionet website.

What I Got

So I downloaded few ECG signals and plotted them in Matlab. Below shows one plot (I lead) I got.

As you can see, this signal is not in the Amplitude 0 line. But all the ECG recordings I have seen and saw in research papers are horizontal. Like in below image.

My Question

I'm new to ECG signal processing. So I'm curious to know,

1)Whether datasets having this kind of abnormal behavior of ECG are wrong?

2)Is there any Matlab function which can convert such data to view like in 2nd image I attached?

UPDATE

To download the dataset,

https://www.physionet.org/cgi-bin/atm/ATM

Go to the above link and select "PTB Diagnosed ECG database" as the Database and select a patient record (E.g: s0010_re).

I have tried the BEADS technique we have developed for the separation of

sparse signals, possibly positive, with sparse derivatives,
a trend (low frequency),
random noise.

I have played with one signal for about two minutes, and here is one outcome:

It is called BEADS, for  "Baseline Estimation And Denoising with Sparsity". The Matlab BEADS code is at MatlabCentral, the paper is Chromatogram baseline estimation and denoising using sparsity (BEADS). It was initially meant for chromatographic signals, but is also used for EEG/ECG and several other spectral signal, see the BEADS page .

One nice feature is that you directly see the noise (and how uncorrelated it can be) and the shape of the trend, then tweak the parameters. Mine (probably not optimal) are:

addpath('C:\Users\duvall\Documents\MATLAB\toolbox\2014_BEADS_Baseline\');
% load('s0010_rem.mat')
 data = val(1,:)';
% Filter parameters
fc = 0.004;     % fc : cut-off frequency (cycles/sample)
d = 1;          % d : filter order parameter (d = 1 or 2)
% Positivity bias (peaks are mostly symmetric)
r = 1;          % r : asymmetry parameter
% Regularization parameters
amp = 0.8;
lam0 = 0.5*amp;
lam1 = 5*amp;
lam2 = 4*amp;

[x1, f1, cost] = beads(data, d, fc, r, lam0, lam1, lam2);
figure(2);
subplot(4,1,1)
h=plot([data]);axis tight;grid on;set(h,'LineWidth',2);
subplot(4,1,2)
h=plot([x1]);axis tight;grid on;set(h,'LineWidth',2);
subplot(4,1,3)
h=plot([f1]);axis tight;grid on;set(h,'LineWidth',2);
subplot(4,1,4)
h=plot([data-x1-f1]);axis tight;grid on;set(h,'LineWidth',2);

Basically, it took me 10-times more  to download and load the signal than to tune the parameters (90 seconds approximately), but I know the code :)


Adaptive filtering: Optimum filter length and delay I'm trying to find the optimum filter length for an Adaptive Filtering, using RLS Algorithm.

I'm using this design:

So the "error" signal is the signal without noise (and that's the signal that I want).

If I have $e(n) = d(n)-y(n)$ but $d(n)$ is my desired signal, I need that $e(n) \rightarrow 0$ so I find the optimum filter length (and the delay) using the MSE criterion, but now I have the signal that I want as the error, so I don't know how to find the optimum filter length, because I have NO idea what MSE do I have to get at the output!

Can anyone tell me what should I do?

Thanks!

In order to be able to choose an optimal value for the delay $\Delta$ it's important to understand how the system works. The purpose of the delay is to decorrelate the desired signal $s(n)$ and the signal component $s(n-\Delta)$ at the input of the adaptive filter. This means that $\Delta$ must be chosen such that the autocorrelation $R_{ss}(k)$ of $s(n)$ is (close to) zero for lags greater than $\Delta$:

$$R_{ss}(k)\approx 0,\qquad |k|>\Delta$$

However, we cannot choose $\Delta$ arbitrarily large because the delayed interference at the input of the filter must be correlated with the interference added to the signal, i.e., the autocorrelation $R_{rr}(k)$ of the interference must still be significant at a lag of $\Delta$, otherwise the adaptive filter cannot predict the interference. If we can assume that $r(n)$ is narrow-band compared to $s(n)$, it's always possible to find an appropriate value for $\Delta$.

With an appropriate value for $\Delta$, the adaptive filter will try to predict the interference, i.e., it will try to undo the effect of the delay in the frequency band where the interference has significant frequency components. So the output of the filter will approximate $r(n)$: $y(n)\approx r(n)$. Consequently, the error signal will approximate the desired signal: $e(n)\approx s(n)$.

After having chosen a value for $\Delta$ based on the autocorrelation of $s(n)$, the filter length must be chosen by trial and error. A long filter will give a better suppression at the cost of slower convergence.


Is it possible to detect or localize movement of objects using LPWA (or LPWAN) such as LoRa, DASH7, NB-FI? Assume that there is a moving human in the room without any device (or receiver). Is it possible to detect or localize movement of objects using LPWA (or LPWAN) such as LoRa, DASH7, NB-FI by analyzing the reflected signal?

Yes, the Doppler effect applies to any electromagnetic wave. No matter whether that is an actual Radar pulse, or a microwave oven falling from an office tower, or a passenger aircraft's transponder, or an LPWAN device.

But, this possibility is especially limited for LPWAN devices, as I'll explain below.

Low Doppler shift due to low carrier frequency

The higher the frequency, the more Doppler, and for some reason you picked particularly low-frequency standards, all operating 1 GHz, making the Doppler shift you get smaller, compared to other uncertainties like local oscillator frequency errors of the typically cheap devices.

Need for Frequency estimation

Short pulses make estimation hard

There's a basic mathematical truth that for a frequency estimate to be fine, you need a long observation, and conversely, for a time estimate to be fine, you need a large bandwidth.

Now, these standards use short transmissions, bursts, because the batteries won't allow them anything else. So that severely limits the time, or energy, you get to estimate frequency. Since you probably don't know in which general direction your LPWAN devices are moving, you'd want an unbiased estimator. Armed with a relatively simple SNR number, you can plug that in to the Cramer-Rao Lower Bound to see how good your frequency estimator of choice will get. And from that it's just multiplication with the speed of light to get the speed estimate variance lower bound. You probably won't like the results.

Cheap devices with no ongoing network registration have bad frequency accuracy

For Doppler estimation to work, you need to be able to estimate the frequency shift of a device. To achieve that, either the device needs to transmit at a very precisely known frequency, by either having a very good oscillator itself, or by "learning" the right frequency from a base station.

Both principles are pretty much impossible for LPWAN devices; they need to be cheap and low-power, so no high-quality oscillator that could run continuously. Also, no "waking up, synching to a base station, correcting the frequency offset one has by complicated algorithms, sending three bytes, going back to sleep", so also the second option is out.

Location estimates hampered by low bandwidth

So, in absence of good transmit frequency knowledge, you could simply observe the same transmitter from multiple points and estimate the speed of movement and the direction through getting enough equations to compensate for the unknowns, right?

Yes. In essence, you'd be building the "inverse" of GPS. Problem: GPS is especially designed to make getting the pseudorange and frequency estimates possible. It sends a relatively wideband, pseudo-white sequence, to make it easy to estimate a sharp crosscorrelation peak.

Your systems are designed to be narrowband and pretty much invariant to static frequency shifts, and hence (see point "short pulses..." above) do not allow for accurate time or frequency recovery.

You basically always run into the fact that these are the worst signals for frequency and/or location detection.

General remarks

I've never heard of DASH7 and NB-FI, but realize the first is just an abused RFID standard, and the second exists as invention of a single company, like there's dozens of companies who just wanted to grab some cash on the IoT goldrush, so honestly, probably not the relevant standards in the future.


Why do we need multiple layers in each octave and multiple octaves in SIFT? I skimmed through the SIFT paper. I understand that there are multiple octaves, which are composed of multiple layers. The layer $k$ of an octave (btw, where does this name come from?) corresponds to the image filtered with a Gaussian kernel with a standard deviation of $k\sigma$, so layers higher in the octave correspond to more blurred versions of the original image. From these blurred images, you calculate the difference of Gaussians (DoG), from which you can compute the pixels that are locally higher (compare to a neighborhood around that pixel). The pixels that are locally higher than its neighbours are the potential keypoints. Once this is done, you can downsample the Gaussian-filtered images in the first octave to produce another octave, then you compute the DoGs again. This process can be repeated (for several iterations).

Intuitively, we apply the Gaussian kernel with a progressively higher standard deviation so that we can find keypoints at different scales (or resolutions), but why do we need to further downsample the Gaussian-filtered images in the first octave and repeat the process? Intuitively, this could still help to find keypoints at different resolutions (or scales), but I can't fully understand why would this approach be correct. So, why do we need multiple layers per octave and multiple octaves in SIFT?

In music theory, an octave is an interval in frequency, from a frequency $f$ to frequency $2f$. For example "an octave higher" means "twice the frequency". Expressed as wavelength inversely proportional to frequency, $\lambda \propto \frac{1}{f}$, an octave would be the interval from a $\lambda$ to $\frac{1}{2}\lambda$. In the SIFT paper's notation, $\sigma$ is the scale, which is analogous to wavelength as it is proportional to the width of the Gaussian kernel.
Downsampling is just to reduce the image size so that there is less computation to do. From the paper (emphasis mine):

Once a complete octave has been processed, we resample the Gaussian
image that has twice the initial value of $σ$ (it will be 2 images from
the top of the stack) by taking every second pixel in each row and
column. The accuracy of sampling relative to $σ$ is no different than
for the start of the previous octave, while computation is greatly
reduced.

What enables to do this? A Gaussian function has another Gaussian function as its Fourier transform, and that decays faster than can be counteracted by any differential operator such as the Laplace operator. A finite-amplitude low-pass signal filtered by a Gaussian or a Laplacian of Gaussian or its approximation of scale $\sigma$ has in an approximate sense a band limit that is proportional to $\frac{1}{\sigma}$. By increasing the constant of proportionality, the mean square power of the part of the signal that is above the band limit can be made arbitrarily small. This allows to resample the filtered signal with a sampling frequency proportional to $\frac{1}{\sigma}$ while keeping the aliasing error arbitrarily small. A convenient approach is to not sample every signal filtered at scale $\sigma$ by sampling frequency $\frac{\sigma_0}{\sigma}$ with $\sigma_0$ some reference scale that is sampled at sampling frequency $1$, but instead by sampling frequency $2^{\operatorname{ceiling}\left(\log_2(\frac{\sigma_0}{\sigma})\right)}$:
Figure 1. Sampling frequency (red) is the same for all frequencies within an octave. The drawback is that for some scales the sampling frequency is higher than the actual requirement (blue) for that scale. $\sigma$ is the scale and $\sigma_0$ is a reference scale.
This bunches together values of $\sigma$ into a number of octaves such that within each octave the sampling frequency is a power of 2, ensuring that the location of each spatial sample conveniently coincides with a spatial sample one octave higher. The drawback is that the sampling frequency for all scales in an an octave must be chosen based on the lowest scale included in the octave, to keep the aliasing noise within acceptable limits for that scale. The higher scales in the octave will typically have less aliasing.


Argument of Dirac delta function is quadratic with complex roots For example, we write $$\delta(at + b) = \frac{1}{|a|}\delta\left(t+\frac{b}{a}\right)$$
In a similar way, $$\delta(t^2 + t + 1)=\underline{\hspace{1cm}}$$

Thanks for the help.

Purna.

I shall provide more details if I am correct. IMO if there is not root on the domain of integration, and here I suppose that $t\in \mathbb{R}$, then the argument never vanishes. Then, in an engineer fashion, one could say (but $t\mapsto \delta(t)$ is not a function):

$$\delta(t^2+t+1) \equiv 0\,.$$

In more precise words, I would consider that for any suitable $f(t)$, $t\in \mathbb{R}$:

$$\int_{-\infty}^{\infty}f(t)\delta(t^2+t+1)dt = 0\,.$$

Related links:

Dirac $δ(g(x))$ with complex roots $x_i $
Regarding $x^2−a^2$ inside the argument of Dirac delta


Convolving With half-Hanning The goal is to get a 200ms decaying delay of an audio signal while preserving the sharpness of attacks to mimick human perception of sound. The paper I'm following convolves a 200ms half-Hanning window with each frequency band to simulate it.

I'm using numpy.hanning() since there's no half-Hanning and have tried either setting the first half of the hanning window to 0, the first half to 1, or chopping it off and only using the second half. I'm also normalizing the window to keep scales the same.

Since I don't really know, which one of these would be the "correct" half-Hanning?

The chopping off version, which seems the most likely to me, works but is offshifted 100 ms to the left. Should I just shift it 100ms to the right and call it good?

dataMasked = np.zeros(data.shape)
for i in xrange(freqs.size): # Loop through each frequency band
    winSize = round(0.4/(length/bins.size)) #Calc win size
    hann = np.hanning(winSize)
    halfPoint = np.argmax(hann)
    hann = hann[halfPoint:] #Chop at half
    hann = normHann(hann)
    dataMasked[i,:] = np.convolve(data[i,:],hann,mode='same')

Here's the plot of what I'm getting (blue is raw loudness curve, while green is adjusted).

The paper I'm following shows this:

Thanks!

Using the second half only (chopping) is the correct way. Note that setting the first half to zero is actually the same thing, up to a constant delay. The 100ms delay you observe might be due to the use of the "same" flag for the numpy convolution. It seems to me that the right truncation for this application would be to use numpy.convolve(signal, window, 'full')[:len(signal)] - without the kind of centering that 'same' does.

As you have noticed, the peaks in the smoothed envelope function do not exactly coincide with the peaks in the original envelope and there is no magic trick to make them coincide (the delay will depend on the sharpness of the attack) - so the resulting curve is not ideal for very accurate beat detection and chopping.


How to calculate transmission capacity? My textbook doesn't explain how they got the answer to this question:

"If a transmission channel uses an analog carrier wave that is modulated by eight amplitudes and sixteen frequencies per bit time, then what is the transmission capacity of that channel in bits per bit time?"

I don't understand how they got the answer and I can't find any good sources online to explain it.

The question as stated doesn't make sense. In particular, "sixteen frequencies per bit time" doesn't mean anything. So, I'll try to guess at what the question actually is, and try to point you in the right direction.

Let's break the question in two parts. First, we're told that the carrier can have eight different amplitudes. This means we can encode three bits in the carrier amplitude. If we label amplitudes $L_0,\ldots,L_7$, we can encode the bit sequence $000$ as amplitude $L_0$, $001$ as amplitude $L_1$, etc. This means that this carrier can transmit bits 3 times faster than a carrier with only two allowed amplitudes.

The second part of the question seems to indicate that you actually have sixteen carriers, all operating simultaneously. (and that they don't interfere with each other). Each carrier can have eight different amplitudes.

If you take $R$ as the baseline rate obtained by a single carrier with two amplitudes, you should now be able to calculate the rate of 16 carriers with 8 amplitudes each, as a multiple of $R$.

Edit: I think there is a second possible interpretation of the question: this could be some sort of combination ASK-FSK, where there are 16 possible frequencies (FSK) and each could have 8 possible amplitudes (ASK). Each frequency would carry 4 bits and each amplitude a further 3 bits, for a total of seven bits per symbol time.


Canny non-max suppression intro I have rich experience signal-processing, but haven't touched image processing since my happy university days... so:

Working on a side project, I need to enhance ridges (continuous lines) in a 2D grey-scale image. Using Matlab, I'm at a point where Canny edge detection does a great job, but since many lines are thick, I get 2-edges for each ridge. Some Questions:

Is Canny a tool fit for the task?
If so, How do I adapt the Canny algorithm to enhance ridges, instead of edges?
Where can I find a good intro to Canny edge detection (googled it, nothing really helpful found)?

Canny's original paper addressed the issue of ridge edges and roof edges in addition to step edges, but I'm unaware of any implementation of "Canny" edge detection that detects anything other than step edges.  Canny's paper is a pretty easy read, but, unfortunately, paywalled.  J. Canny, "A computational approach to edge detection", IEEE Trans. Pattern Analysis and Machine Intelligence, 8:679-714, 1986.

Canny's paper is really about how to use calculus of variations to derive a 1-dimensional filter that has optimal response to each type of edge.  It turns out that the ridge-edge filter is not scale invariant (you have to know how wide your ridges are).  Here is the picture from Canny's paper of the optimal ridge and roof detection operators (page 683, © IEEE, 1986, this is "fair use"):

What is typically called "Canny edge detection" is as follows:

Perform an approximation of a derivative of Gaussian function in X and Y.  (This is usually done with a Gaussian blur followed by the Sobel operator.)  (The derivative of Guassian is actually not the optimal step-edge detection operator, but the optimal step-edge detection operator "kinda-looks-like" a derivative of Gaussian.)
Convert from $\langle\Delta x, \Delta y\rangle$ form at each pixel to $\langle$magnitude,"angle"$\rangle$.  (Where "angle" is a 2-bit approximation that tells you whether the edge is closer to horizontal, vertical or one of the two 45-degree diagonals.)
"Thinning": For each pixel look at its two neighbors in the "angle" direction and keep the pixel only if it is the maximum of itself and its two neighbors.  (Otherwise set the magnitude to 0.)
"Thresholding": For each pixel with magnitude less than some "weak threshold" suppress the pixel to 0.  For each pixel with magnitude greater than some "strong threshold" set the pixel to whatever value means "definitely an edge."  For pixels between the two thresholds set the pixel to the value "maybe an edge".
Hysteresis: For each "maybe an edge" pixel change it to "definitely an edge" if any of its 8 neighbors in the horizontal, vertical or diagonal directions is "definitely an edge."

Theoretically, you should be able to turn a Canny edge detector into a ridge-edge detector by changing just step 1, but I've never tried it.


Impulse response of a time scaling system Assume a bandlimited signal $X(t)$. Given that the output for this signal is $X(t/2)$, what will be the impulse response $h(t)$ of such a system?

\begin{array}{l}
X( \omega ) \ =\ \int ^{\infty }_{-\infty } x( t) e ^{-j\omega t} dt  \\
Y( \omega ) \ =2X( 2\omega )\\
\end{array}

As the system won't be LTI the output won't be the convolution of the input and the system response. How do I proceed with finding the response of such a system? Any guidance in the direction of the solution will be very helpful.

Edit: Thanks for all the answers, i had reframed this question which was on a similar line. The question goes like this ->
Given an input pulse of duration T, the system scales it to a pulse of 2T, design such a system? It would be of great help if someone could point out the right direction to proceed.

The input-output relation of such a system can be written as

$$y(t)=\int_{-\infty}^{\infty}\delta(t/2-\tau)x(\tau)d\tau\tag{1}$$

Note that the system is linear but time-varying, and such systems can generally be described by a two-dimensional impulse response $h(t,\tau)$:

$$y(t)=\int_{-\infty}^{\infty}h(t,\tau)x(\tau)d\tau\tag{2}$$

In the given example we simply have

$$h(t,\tau)=\delta(t/2-\tau)\tag{3}$$

For linear time-invariant (LTI) systems, $h(t,\tau)$ only depends on the difference $t-\tau$.


Filter before or after multiplication of two signals? I have two signals and want to multiply them, e.g. voltage and current are multiplied to get power. The result (e.g. power) shall be a filtered result, e.g. to see the average.

Is it better to filter (digitally) each signal separately before the multiplication or  to filter after the multiplication?

Intuitively, I would have said that filtering before the multiplication is better for three reasons:

Less numerical problems
Noise of both signals is added by the multiplication, thus filtering before multiplication improves SNR
Multiplication generates modulation products and if you filter out a frequency you do not want to see anyway it also does not cause modulation products.

Am I on the right track here or is there something I am overlooking?

Related: SNR After Multiplying Two Noisy Signals

Edit 1: To answer the question of typical signal spectra: Think of a laboratory grade measurement device like a multimeter. It acquires two signals which can be any signal shape (in time domain) but are band-limited by the input circuitry to avoid anti-aliasing. Noise is typically white noise with a bit of 1/f noise added. SNR is usually 10's of dB.

A personal rule: in general, it can be better to perform non-linear operations before linear ones. One reason behind that is that a lot of practical concerns are related to outliers or suspect behavior, which can  easily be smoothed out (and become indistinct from other signals) by  linear filters.

Let me reformulate. If $f_i$ denote filters, and $s_i$ signals, should one do $f_0 \ast(s_1 .s_2)$ or $(f_1 \ast s_1)(f_2 \ast s_2)$? Without a better knowledge on the spectra of the signal, the nature of the noise and the numerically objective, I have no definite answer. The second one is apparently more flexible, become you can play on two filters.

If the result is some kind of average, I would filer the product, as  averaging is filtering, and it commutes with filtering. Thus, you have more chance to preserve the average after filtering the product, if the noise is low. Since two signals like $+1,-1,+1,-1\ldots$. Any even-length avering will turn them into $0$, while their product is flat with average 1.

I cannot see why you will have less numerical problems.
Noise is added "logarithmically", but the same happens to signals. But yes, modulation can be troublesome, but two cosine with zero mean will be multiplied into a non-zero mean signal, and this is what you want to measure.

In some cases though, to preserve contrast in image, homomorphic filtering convert a product into a sum to filter separately the two components (cf. "The second one is apparently more flexible"). Finally, you can use non-linear filters on the product as well, in case of complicated noises.  A related question on the ordre of operations can be found in Is convolution distributive over multiplication?


Design of Parks-McClellan Lowpass Filter Can we design Parks-McClellan Lowpass Filter with only sampling frequency of 12KHz and cutoff frequency of 4KHz and order of 30? I went through Matlab help, where we need parameter like passband cutoff frquency, stopband cutoff frequency, passband ripple e.t.c.,

The minimal specification for a FIR lowpass filter designed by the Parks McClellan algorithm is the passband edge, the stopband edge, and the filter order. You don't need to specify passband and stopband ripples; they are the result of the design process. What you can do is specify passband and stopband error weights, which let you specify the relation between passband and stopband approximation errors (a larger weight reduces the error).


Is $y(t) = \cos(t) + x(t)$ a time-invariant system? Is $y(t) = \cos(t) + x(t)$ a time-invariant system?

$y(t-k) = \cos(t-k) + x(t-k)$

But it isn't equal to $\cos(t) + x(t-k)$

So, would it be time-invariant?

You've answered the question yourself already. So, does the response to a shifted input equal the shifted response to the original (non-shifted) input? If so, then the system is time-invariant. Otherwise, it's a time-varying system.


How to remove 50 Hz frequency noise from a PCM file? Are you kind to tell me the odds or the means by which should I grab a clean signal with the $~50 \mathrm{Hz}$ hum removed using just programming?

I managed to make a program in Visual-C that composes a sinusoid and overlaps it over a signal in a file. The file is recorded using a microphone. At recording time using a PC I get the frequency hum of the AC, with the phase with which it overlaps the signal. I see that the noise's frequency is not influenced and remains constant over time.

You should just use a notch filter.

Below shows two sinusoids added together (top) and the result of filtering the lower frequency one (bottom).

The filter is just:
$$
H(z) = \frac{1 - 2 \cos(2\pi \omega_{b}) z^{-1} + z^{-2}}
{1 - 2 \alpha \cos(2\pi \omega_{b}) z^{-1} + \alpha^2z^{-2}}
$$
where $\omega_b$ is the bad frequency and $\alpha<1.0$ is related to the notch width / depth.

R Code Below

#29936

fs <- 8000
fbad <- 50
phibad <- 2*pi*runif(1,0,1)
fgood <- 101.98340234
phigood <- 2*pi*runif(1,0,1)

T <- 1000
t <- seq(1,T)

x_noisy <-  sin(2*pi*fgood/fs*t + phigood) + sin(2*pi*fbad/fs*t + phibad)

num <- c(1, -2*cos(2*pi*fbad/fs), 1)
alpha <- 0.99
den <- c(1, -2*alpha*cos(2*pi*fbad/fs), alpha*alpha)

x_filtered <- filter(num, den, x_noisy)

par(mfrow = c(2,1))
plot(t,x_noisy, type="l")

plot(t,x_filtered, type="l")


Wiener Filters - Obtaining / Estimating Image and Noise Spectrum We all know the famous Wiener Filter,

$$G(\omega_1 , \omega_2)=\dfrac{H^*(\omega_1 , \omega_2)}{|H(\omega_1 , \omega_2)|^2 + \dfrac{N(\omega_1 , \omega_2)}{S(\omega_1 , \omega_2)}}$$

How to we use it practically? Meaning, given $y = s * h + n$ where $s$ is the original image, $h$ is a LSI operator, and $n$ is some noise (independent of $s$), how to calulate (or estimate) $N(\omega_1 , \omega_2)$ and $S(\omega_1 , \omega_2)$ in order to reconstruct $s$ from $y$?

You are speaking about so called Wiener deconvolution, because Wiener filter is set by equation $ w = R^{-1}r $, where:

$ R = E[y(n)y(n)^H] $ - autocorrelation matrix of input signal
$ r = E[y(n)s(n)^*] $ - crosscorrelation vector of input signal and original signal.

So the practical solution when statistics $R$ and $r$ are unknown is the Method of Least Squares, where only data matrix $F$ and part of original signal $s$ are needed.

In your case you should have some a priori knowledge about formula entities. Try:

  https://www.umanitoba.ca/faculties/medicine/radiology/Dick_Gordon_papers/Dhawan%20et%CAal.%20%281985%29.pdf

or:

  http://www.eos.ubc.ca/research/cdsst/members/6_96.pdf

It's a pity I only concern this subject in radio applications so there is nothing more to advise.

Hope this helps!


Can I take piecewise FFT of a larger input sample? I have a N point sequence x. Let X be the fft(x). (Assume x is arbitrary)

where,

x=[x1 x2 x3 x4 x5 x6 x7 x8]

fft(x)=[X1 X2 X3 X4 X5 X6 X7 X8]

How do I get the result of a 8 point DFT by doing cascading 4 point DFT and 2 point DFT?

Suppose I reshape the signal 'x' (of dimensions '1 x 8') to 'x1' (of dimensions '2 x 4')

x1=[x1 x2 x3 x4;
   x5 x6 x7 x8]

fft(x1 (first row)) i.e. fft(x1 x2 x3 x4) say =  [Y1 Y2 Y3 Y4]
fft(x1 (second row)) i.e. fft(x4 x5 x6 x7)  say = [Y5 Y6 Y7 Y8]

Then again taking fft column wise

fft(Y1 Y5)  say = [Z1 Z2]
fft(Y2 Y6)  say = [Z3 Z4]
fft(Y3 Y7)  say = [Z5 Z6]
fft(Y4 Y8)  say = [Z7 Z8]

However this Z will not be equal to X. How do I proceed to get X by doing these smaller point FFTs?

Any help would be great! Thanks

You need the Cooley-Tukey algorithm, which can be used to express the DFT of size $N=N_1N_2$ by DFTs of sizes $N_1$ and $N_2$. It works a bit like what you sketched in your question, but you forgot the twiddle factors. This answer explains the details of the algorithm and it also includes a simple Matlab implementation.

EDIT (see comments): If you really want to compute the FFTs of blocks of consecutive data you need to split the computation of even and odd frequency indices (assuming the DFT length $N$ is even):

$$X[2k]=\sum_{n=0}^{N/2-1}\left(x[n]+x[n+N/2]\right)W_{N/2}^{nk}\\
X[2k+1]=\sum_{n=0}^{N/2-1}\left(x[n]-x[n+N/2]\right)W_N^nW_{N/2}^{nk}\tag{1}$$

where $W_N=e^{-j2\pi/N}$. Note that you can of course compute the DFTs of the length $N/2$ sub-blocks separately instead of adding (or subtracting) the blocks before the DFT (as in (1)), but this is computationally less efficient.


Generay way to find Stability of system I am trying to find a general way to find system stability.I have applied these methods and struck in confusion
For eg:
If
$$y(t)= \int_{-\infty}^{t} x(\tau) \sin(4\tau) d\tau$$
then find whether the system is stable or unstable?

Process (1):
If $$y(t) = \int_{-\infty}^{t} x(\tau)\sin(4\tau) d\tau$$
then
$$h(t) =  \int_{-\infty}^{t} \delta(\tau) \sin(4\tau) d\tau$$
$$\implies h(t) =  \int_{-\infty}^{\infty} \delta(\tau)\sin(4\tau)  u(t-\tau) d\tau$$
$$\implies h(t) = [\sin(4\tau)u(t-\tau)]|_{\tau=0}$$
          $$\implies h(t)=0$$.

so system is stable as $$\int_{-\infty}^{\infty} |h(t)| dt= 0 < \infty$$ ;

Process (2):
If $$y(t)= \int_{-\infty}^{t} x(\tau)\sin(4\tau) d\tau$$
then when $x(t)=\sin(4t)$
$$y(t)= \int_{-\infty}^{t} \sin(4\tau)\sin(4\tau) d\tau  = \infty $$
so the system produces unbounded output $y(t)$ for bounded input $x(t)= \sin(4t)$
Hence the system is unstable

Thus actually what will be the system stability and which process is right and which is wrong and how it is wrong?? please explain....

Your first method assumes that the system is entirely characterized by its response to the input $x(t)=\delta(t)$. This would only be the case if the system were linear and time-invariant (LTI). In that case, the input-output relation would have the form

$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau\tag{1}$$

where $h(t)$ is the system's impulse response. However, the given input-output relation is not of the form $(1)$, and hence the system is not LTI. Consequently, the system's response to $x(t)=\delta(t)$ does not say much about the system's behavior in general. That's why your first method of determining the stability is not valid.

The second method is valid. You found a bounded input signal for which the output is unbounded, and in this way you've shown that the system is not BIBO-stable.


Wavelet image denoising: dual-tree versus double-density What is the main difference between the dual-tree DWT and the double density DWT  image denoising techniques?

How does the combination of them improve the quality of the denoised image?

The formalism of the continuous wavelet transform is relatively flexible. To make a practical tool out of it, you ought to discretize it, and here come the pain. It is quite easy to discretize it in a very redundant way (a discrete wavelet frame), but the price to pay when analyzing high-dimensional signals (2D images, 3D volumes) can quickly become too high: on the order or $J^d$, with $J$ the number of wavelet levels and $d$ the dimension.

So people have tried to obtain a critically sampled, discrete wavelet, without redundancy. This is the discrete wavelet transform (DWT). Nothing comes free. The price to pay is a narrow choice in the wavelets, and a processing sensitive to shifts. Ivan Selesnick, among others, worked on several intermediate designs: by allowing a bearable redundancy (typically $2\times$ or a little more in each dimension), one can obtain results close to, or even above the quality of fully redundant designs.

He worked especially:

on dual-tree wavelets and
on the double-density wavelet,

both  with a redundancy of $2$ in 1D, whatever the level. The first one compensates the poorer DWT choice by combining it with its Hilbert pair, allowing lower shift variance. The second allowed more degrees of freedom in the design of wavelet filters, hence a richer choice. Both have their merits, that may depend on the application. A little redundancy in wavelets often helps.

And both can be further combined, as in  The Double-Density Dual-Tree Discrete Wavelet Transform, by the same author in 2004. The gain comes from the double advantage  of having more design freedom in the wavelets (allowing for instance better cutoff frequency and symmetry) and the Hilbert pair for shift-invariance.

This is confirmed at Double-Density Complex Wavelet Transform:

  Both the double-density DWT and the dual-tree DWT have their own
  distinct characteristics and advantages, and as such, it was only
  natural to combine the two into one transform called the
  double-density complex (or double-density dual-tree) DWT.


Does Discrete Wavelet Consider Sampling Frequency or Cutoff Frequency of Bandpass Filter? I have EEG signal with a sampling frequency of 128 Hz, and filtered from 1-30 Hz using a bandpass filter. Now I want to know while decomposing the signal; does DWT decomposes it based on sampling frequency or cutoff frequency? Detailed Coefficients one should be 7.5-30 (if consider cutoff frequency), or 32-64 (if consider sampling frequency).
If we consider the cutoff frequency, we now that upper limit is a bit higher than 30, so how can we correctly justify it?
I am using python pywavelets and want to compute the frequency range of decomposed level.

THe DWT does not really know about the actual sampling. It do cares about the relative frequency span. And it has not knowledge about pre-processing, like linear filtering.
So in your case, the detail coefficients of the first level of wavelet coefficients correspond,  more or less, to the [128/2-128]/2 frequency range.
Because of the initial [0-30] low-pass filter, the above [32-64] first wavelet subband range is likely to contain little useful information.


What will be the distribution of a given signal Suppose I have a signal $$
\mathbf{s} \sim C \mathcal{N}\left(\mathbf{0},  \sigma^{2} \mathbf{I}_{N_{}}\right)
$$
The signal going through a channel :
$$
{\mathbf{r}}_{}=\mathbf{s}+{\delta \mathbf{n}}
$$
where :
$$
\mathbf{n} \sim C \mathcal{N}\left(\mathbf{0},  \sigma^{2} \mathbf{I}_{N_{}}\right)
$$
$$
\mathbf{\delta} \sim  \mathcal U\left(-1, 1\right)
$$
What will be the distribution of ${\mathbf{r}}_{}$. I am having difficulty because of the second term where the uniform distribution scalar multiplied with a Gaussian vector.

The way to proceed is to use the definition of Marginal distribution, and also Law of the unconscious statistician,
$$p_\mathbf{R}(\mathbf{r}) = \mathbb{E}_\delta\big\{p_{\mathbf{R}\mid \delta}(\mathbf{r}\mid\delta)\big\} = \int_\delta p_{\mathbf{R}\mid \delta}\left(\mathbf{r}\mid\delta\right) p_\delta(\delta) \mathrm{d} \delta = \int_{-1}^{+1} p_{\mathbf{R}\mid \delta}\left(\mathbf{r}\mid\delta\right) \frac{1}{2} \mathrm{d} \delta \tag{1}$$
Note that for a given $\delta$, the conditional $\mathbf{r}$ is $C\mathcal{N}\left(\mathbf{0}, \sigma_s^{2} \mathbf{I}_{N} + \delta^2 \sigma_n^{2} \mathbf{I}_{N}\right)$.
The remaining work is to inject the conditional pdf to the uttermost right hand side of $(1)$ and calculate the integral.
An elegant solution would be using the nice properties of Gaussian random variables: sum of (independent) Gaussian is also Gaussian. In that way, one can push the integral to the covariance matrix part and come up with the result:
$$C\mathcal{N}\left(\mathbf{0}, \sigma_s^{2} \mathbf{I}_{N} + \frac{1}{3} \sigma_n^{2} \mathbf{I}_{N}\right)$$
Bonus. Use this script to check the result of the scalar case.
sigmas = 2; sigman = 5; N = 1e5;
s = sigmas/sqrt(2)*complex(randn(N,1),randn(N,1));
n = sigman/sqrt(2)*complex(randn(N,1),randn(N,1));
delta = 1-2*rand(N,1);
r = s + delta.*n;
histogram(real(r),'Normalization','probability')
xx = linspace(-10,10,40); yy = normpdf(xx,0,sqrt((sigmas^2+sigman^2/3)/2));
hold on;
plot(xx,yy,'ro');


How is noise added in an image during image acquisition or image transmission? Recently I read about image noise. I read about noise characteristics, but I am unable to understand the following two points:

The sources of noise in digital images arise during image acquisition
(digitization) and transmission.

Imaging sensors can be affected by ambient conditions
Interference can be added to an image during transmission

Imaging sensors are just electronic devices that receive photons (light) and turn them in to electrons (electric signal). As such, they are affected by ambient conditions such as temperature.
For example, you will have some noise coming from the electronic agitation (false currents due to freely moving electrons, without incoming light) that grows with the sensor temperature.

You do also have non-linear responses in the sensor that cause other noise patterns. A pixel needs to receive enough light before providing a reliable measurement. Before that threshold, it will be affected by Poisson noise (a noise whose distribution depends on the value of the pixel). On the other side of the per-pixel energy, it can get saturated (too much light), and its energy will "spread" onto neighboring pixels, thus modifying their value.

Finally, since you have a digital device, then you have an analog-to-digital conversion step that is likely to introduce further quantization noise, i.e., the recorded values will not be the "real" ones but something close up to the measurement precision.

Outside your imaging device, there is light transmission in the atmosphere (turbulence is very likely to be a disturbance cause here, but it's usually not considered as noise), and image transmission over communication channels where packet losses (and thus pixel loss!) can happen.


data rate with a FSK modulation With a BPSK modulation the data rate is $D=\frac{1}{T_s}$ where $T_s$ is the symbol time. But for a BFSK modulation, how do we find the data rate?

In general, the data rate in bits per second is the symbol rate times the number of bits per symbol; independent of the modulation used.

The B in both BPSK and BFSK stands for binary, so there is 1 bit per symbol, or 2 possible symbols (2 possible phases or frequencies depending on the case).

So the data rate in both cases is the same.


Detection of "constant" signal I have the positions of a computer mouse, sampled at a high frequency, with a good amount of noise. When the user stops moving the mouse, instead of going to zero the velocity reported by the mouse is non-zero and the cursor drifts.

I want to estimate the drift velocity. The attached image is a plot of the x and y velocity components. The user has is not moving the curs between datapoints ~600 and ~660.

Are there standard techniques for identifying a constant subinterval of a signal?

Here are some things I have tried:

Using a high pass filter to estimate the standard deviation of the noise. Then I used that as a comparison of nearest neighbor velocity differences. This does not catch slowly increasing intervals. (Black intervals in second figure)
Looked at instantaneous min/max over the current interval and filtered out those intervals where (max-min) is too high. (Red intervals in second figure).

As you can see from the second image, I can select out some reasonable intervals of interest, but I still end up with spurious ones. I could come up with more filters, but I'm now just introducing more ad hoc parameters.
What is a more principled way to do this?

Perhaps transferable to your problem, one method from Olyha B., Rutledge J., "TrackPoint System Version 4.0 Engineering Specification", IBM, 1999:

  2.2.3 DRIFT CORRECTION

  Due to the significant temperature sensitivity of the force sensors, it is necessary to recalibrate the zero force
  origin of each axis on a periodic basis. In order to properly
  recalibrate, it is necessary to determine when a ‘hands off’ condition
  exists (i.e. when no force is being applied). This is done by
  monitoring the X and Y axis force samples, and when there is minimal
  change in the readings (no more than ±1 count) for a specified length
  of time (3 seconds by default) it is assumed that the hands off
  condition exists. At this time the X, Y, and Z axis origins are reset
  to the smoothed current value of the X, Y, and Z axis force samples
  respectively. It should be noted that noise introduced by any source
  other than applied sensor forces will severely impact the
  effectiveness of this technique.

In my experience as a long-time ThinkPad user, this technique cannot always differentiate between the user trying to move the mouse cursor at constant speed and real sensor drift. Then again this is not usually dangerous and it is easy to manually trigger recalibration by hands off. Perhaps there should be a more advanced analysis algorithm running that can identify velocity variation in a manual attempt to move the mouse cursor at constant speed.


Frequency filtering - Audio Sound are digitalized in computers as an archive with its intensity amplitude at the given time. I want to ask a very simple question about frequency filtering. How it's done, if there's no frequency information stored in the digital audio? It's all intensity. I cannot understand.

An set of digitized intensity samples that varies in time produces a waveform (just graph the data).  If the waveform shows any periodicity, that will be an obvious frequency.  Even if it doesn't, Fourier theory says any (non-pathological) waveform can be decomposed into a set (perhaps infinite) of sinewaves.  And sinewaves are describe as having a frequency.  Thus frequency information is implicit.

And this frequency information can be changed by certain filtering processes or algorithms many different ways.

Even if the digitized samples don't vary at all (a set of constants), that is the same as DC, or a frequency of zero (0 Hz).


How can I calculate the height at which the mobile phone is located? I need to calculate Hm - the height at which the mobile phone is located.
For calculating I can use:

Hb - the BTS height, m
L - the rssi-parameter, dB,
f - frequency, MHz.
When I use an Okuruma-Hata model's expression , I get abnormal data.

May be you know any other expressions for calculating the height at which the mobile phone is located?

That model is used to estimate path loss based on distance and frequency (applicable to RF frequencies from 150 to 1500 MHz). Since there are an infinite number of heights that would resolve to the same distance, you cannot use the formula with the transmission from a single BTS to resolve the mobile height.
You could theoretically triangulate this with multiple transmissions from multiple base stations but given the variance in the actual land-based propagation models there would be significant challenges in achieving any sort of reasonable accuracy.
Why not use GPS since nearly every mobile has a GPS receiver now?


Impulse response of a linear time-variant system Is the impulse response of a LINEAR TIME-VARIANT system unique?
I know that it's unique for the L.T.I case!!!
Advance thanks for the HELP....

With the input-output relationship of a linear time-varying system given by

$$y(t)=\int_{-\infty}^{\infty}h(t,\tau)x(t-\tau)d\tau\tag{1}$$

it is straightforward to show that the kernel $h(t,\tau)$ (also referred to as impulse response or input delay spread function), must be unique.

Uniqueness means that there is no other kernel $g(t,\tau)\neq h(t,\tau)$ such that the output $y(t)$ is the same for all possible input signals $x(t)$.

Take the input signal $x(t)=\delta(t-t_0)$. The corresponding output signal is $y(t)=h(t,t-t_0)$, which, by definition, is different from the output signal $y(t)=g(t,t-t_0)$ obtained with the kernel $g(t,\tau)$.


Why are 2nd order edge detectors such as Laplacian more susceptible to noise? I've been told that they are influenced by noise considerably more than that of a first order edge detector such as a Gaussian. But why?

Generally 2nd derivative is more sensitive to noise than 1st derivative. The 2nd derivative is usually accompanied by zero crossing detection, so it works better when grey level transitions are smooth. When the noise level is high, the gray level fluctuates and may result in many zero crossings.

Besides, there is less chance of multiple edge responses with 2nd derivative edge detectors.


Pitch detection to find music notes (from voice) Hi signal processing community,

currently I am searching a topic for my master thesis in computer science and since I am a very good singer I want to do something with music.
I thought about a topic with pitch detection algorithms. But I am very new in this context and I want to inform me first, so I hope you can help me and answer me some of my questions.

I know that I have to work with FFT (Fast Fourier Transformation) and algorithms like AMDF or Yin. Can you give me a hint of a starting point (a paper or website).

I start with a signal which I got from my microphone, but what are the next steps.

I hope you can help me, thank you all.

1) Get some practical experience on what you are studying.  By this I mean is take your recordings of your voice and look at them close up on you computer.  You can use a program like Audacity.

2) Yes, you will have to learn about the Discrete Fourier Transform (DFT, the FFT is a computationally efficient version).  I write blog about it explaining it in my own way.  You can find the link on my profile page.

3) Pick a different topic that is worthy of a Master's Thesis.  There are many techniques for pitch determination, which is a pyscho-acoustic phenomenon vs frequency estimation which is basically a physics problem.  Since you are a singer, you might be interested in how sound waves are attenuated as they travel through air.  You simply have to set up some microphones at different distances from the source and look at their waveforms to see what I am talking about.

4) Depending on your math chops, you may consider coming up with data visualization techniques for signal data.  There is lots of room for creativity here.

Hope this helps.

Ced


What is the proper way to implement a real time spectrogram? I am trying to implement a real time spectrogram and I am not sure that I am implementing it the right way.  First I decided to use Matlab, I was taking a short recording of lets say 0.1 seconds and then perform the FFT on that array of numbers which basically consists of the amplitude of lets say a 11025 sampling rate in that 0.1 seconds.  I was then plotting the result of the FFT after computing the magnitude.  After this process, I would then record another 0.1 seconds and so on.  This was being done in an infinite loop in order to try to achieve a real time result. I know that this is not a spectrogram is just a different representation of this information.  The problem with this is that I am not getting a very good results from this calculation.  For example I was giving a certain tone using a frequency tone generator, and I was seeing a peak at just some instances which shows that somewhere the algorithm is failing.  I do not have to use MATLAB, I could also use other programming languages such as Java so I prefer if I do not utilize a ready made function such as 'spectrogram' which is only present in MATLAB.

So can someone tell me on how a real time FFT of audio should be performed please? So how real time spectrum analyzer achieve that quality and they are real time so if they are taking a short sample of time it must be really short and still manage to get a good representation out of it.

Spectrograms are normally computed using the Short-Time Fourier Transform (STFT), and can be done in real-time. It consists of splitting the signal into overlapped frames, using a window function (typically scaled cosine like Hann) and then computing the FFT for the frame.

The overlap and windowing is key to avoid 'spectral leakage', and sounds like what is missing in your case.


Determine Frequency Resolution samples Can someone tell me if the way I solved this problem is right?

Suppose I have a signal $X_a(t)= 0.5cos(700πt)+0.6cos(720πt)+0.1cos(780πt)$ where $F_s=8000Hz$

I want to determine the minimum samples which must be available in order the frequency resolution of DFT be ok.

I did this:

Firstly calculate each frequency of the signal:
$F_1=350Hz$, $F_2=360Hz$, $F_3=390Hz$
Calculate the differences: $F_2$-$F_1$=$10Hz$ ,$F_3$-$F_2$= $30Hz$ , $F_3$-$F_1$= $40Hz$
$ΔF_{min}$=$10Hz$
$ΔF$ = $ \frac{f_s}{N} $, so $N\geq\frac{8000Hz}{10}$=800 samples

Can someone tell me if my solution is right and explain why we choose the minimum difference

Thanks

It's usually considered difficult to resolve a pair of adjacent frequencies unless there is dip between the two FFT peaks.  Depending on the window (rectangular or otherwise) and how much of a dip you require between the two peaks to consider them resolved (above some floor), and whether or what kind of result interpolation you use, you will need from a bit less than double to triple or more of your suggested 800 FFT result bins.

Of course in a-priori zero noise, your signal of 3 pure sinusoids has 9 free parameters, so you might be able to determine your 3 input frequencies with perhaps as few as 9 or a bit more non-aliased samples.  In any case, much less than 800.

The resolution of a single isolated frequency peak above a low enough noise floor to the nearest 10Hz can likely be done using far fewer than 800 samples by using interpolation methods on the FFT result, depending on the SNR.


Sub-band factor 0.5 in uniform filterbank Is there someone who knows where the factor $\frac{1}{2}$ in front of $U_i(z)$ comes from?

The factor of 1/2 comes from the fact that you are down-sampling. The down-sampling property of the z-transform states that given:

$$ y[n] = x[nL] $$

Then:

$$ Y(z) = \frac{1}L \sum_{r=0}^{L-1}X(z^{1/L}e^{-j2{\pi}r/L})$$

And in your case $L = 2$ so $r$ indexes from 0 to 1 so then:

$$ Y(z) = \frac{1}2 \sum_{r=0}^{1}X(z^{1/2}e^{-j{\pi}r})$$

To be more thorough:

$$ y[n] = u_i[n] = (x * h_i)[nL] $$

So:

$$ Y(z) = U_i(z) = \frac{1}L \sum_{r=0}^{L-1}H_i(z^{1/L}e^{-j2{\pi}r/L})X(z^{1/L}e^{-j2{\pi}r/L})$$

And again $L = 2$ so $r$ indexes from 0 to 1 so then:

$$ Y(z) = U_i(z) = \frac{1}2 \sum_{r=0}^{1}H_i(z^{1/2}e^{-j{\pi}r})X(z^{1/2}e^{-j{\pi}r})$$


Why is the green colour not represented as black in Cr component of YCbCr? I am testing conversion from RGB to YCbCr in Python's pillow imaging module (formerly PIL), and am getting predictable results overall, except for green, as represented in the Cr channel component of the YCbCr version of image.

For blue, I get white (RGB (255,255,255); V in HSV of 100%) in the Cb channel, and black (RGB (0,0,0); V in HSV of 0%) for yellow in the same channel. In Cr, I get white for the red colour, but (RGB (21,21,21); V in HSV of 8%) for green. Shouldn't green be black in Cr?

I thought I would check here before claiming there is a bug in pillow. Perhaps I'm missing something about YCbCr, and the green of it is not the same as the RGB green. But I don't think so, as blue, yellow and red do match RGB.

Found the answer with more experimentation: the opposite colour to red in Cr is not green but cyan, RGB (0,255,255), and it does appear black in the Cr channel. This makes sense, as yellow, RGB (255,255,0), is opposite blue.


Understanding the Inverse DFT (FFT) I was reading this article to understand the fft - JACK SCHAEDLER - A WALKTHROUGH OF THE DFT - STEP BY STEP.
I understand that:

the fft takes a block (256 or 512) of the signal, creates a sinusoidal and cosineidal wave of frequency "f" and multiplies each sample for each sample of the signal to decompose,
add the value of all the samples of the cosine and sine wave product, if it gives greater than 0 it means that the frequency wave "f" is present in the signal
the result of the sum of each product wave gives a complex number

now what I want to know what is the next formula to rebuild the signal, ie an inverse fft

If you know Linear Algebra then the process can be explained using Linear Algebra in a very simple way.

Let's say you have Orthonormal Basis for $ \mathbb{R}^{k} $ given by the Columns of $ A \in \mathbb{R}^{k \times k} $.

Now imagine having a vector $ x \in \mathbb{C}^{k} $.
How would represent the vector $ x $ in the set of coordinate of matrix $ A $, namely using the Orthonormal Basis?

Well, you'd calculate the projection (Inner Product) between any column of $ A $ to the vector $ x $ by $ {x}_{A} = {A}^{H} x $.
How would you recreate $ x $ from the columns of $ A $ using $ {x}_{A} $?
Well, you can see that $ x = A {x}_{A} $ since $ A $ is Unitary Matrix.
Basically you sum the columns of $ A $ using the weights of $ {x}_{A} $.

So, Now just replace $ A $ with the DFT Matrix and you have the whole process of DFT and Inverse DFT.

In practice, the Inverse DFT and the DFT are very similar, multiplication and summation.
The only difference is on one direction (Forward) we multiply by $ {e}^{- \frac{2 \pi j}{N} k n} $ the signal and on the backward we multiply the DFT coefficients by $ {e}^{\frac{2 \pi j}{N} k n} $.


Frequency slows faster when it's decreasing I'm sure there's a simple explanation but I can't discover it. Consider two frequencies, a and b. a is higher than b. I can create a wave with both frequencies. But when I create a wave with frequency decreasing from a to b, the frequency decreases to even slower than frequency b. Consider the following MATLAB code.
timeVals = 0:0.001:1;
sourceFreq = 8;
freqPulse = linspace(1,.5,length(timeVals));
sourceActivity1 = cos(2 * pi * (sourceFreq .* freqPulse .* timeVals));

freqPulse = linspace(.5,.5,length(timeVals));
sourceActivity2 = cos(2 * pi * (sourceFreq .* freqPulse .* timeVals));

plot(sourceActivity1); hold on; plot(sourceActivity2)

Now, compare the two sourceActivities. The end frequency of sourceActivity1 is even slower than 4. And the frequency at the end of the wave is slower than sourceActivity2. This can be determined using visual inspection alone, or with spectrograms.

I think the issue is that you aren't ramping down the frequency from A to B, but the phase.  The phase needs to be integrated from the last previous value.
Plot the phase argument and you'll see what you did wrong.  The phase for source activity 1 should end up parallel but above the phase for source activity 2.
You say that the end frequency is 4 in both cases, but what you really have done is constrict both waveforms to have the same number of cycles.  Since source activity 1 starts faster, it must end slower.


This question is about modelling mechanical response of behavior in terms of signals My dissertation supervisor ask me to code this.

r1=.3;
r2=.5;
xh=[0.00001:.01:.21];
yh=3.2*xh-0.8;
A=2*r2*xh;
B=2*r2*yh;
R=sqrt(A.^2+B.^2);
C=xh.^2+yh.^2+r2^2-r1^2;
z=C./R;
zz=sqrt(1-z.^2);
theta2=atan2(B,A)+atan2(zz,z);
theta1=atan2((yh-r2*sin(theta2)),(xh-r2*cos(theta2)));
plot(xh,180*theta2/pi)
hold on
plot(xh,180*theta1/pi,'r')
axis([0 .2 -170 0])
xlabel('xh [m]')
ylabel('\theta [degree]')
legend('\theta_2','\theta_1')
title 'The change in angles between the elbow and shoulder while lifting a heavy Load' grid

I understand it but now he has ask me to..

" Please do the links velocities and send me the results, please plot velocity of the links against yh, limit the axis of yh to be  ‘axis([-.76 -.16  -15 5])’"

I know that length of r1 or r2 multiplied by theta1 or theta2 should give arc length. Then arc length divided by a time interval will give rad/sec
which i can plot against yh...........I dont know what this axis([-.76 -.16  -15 5])’ means...... BUT I AM VERY NEW TO MATLAB HELP PLEASE.

I dont know what this axis([-.76 -.16 -15 5])’ means

As far as I am concerned axis([-.76 -.16 -15 5]) set minimal and maximum values for x axis (first and second values) and y axis (third and forth values) in plot.

For example:

x = linspace (-10, 10, 100);
y = x .^ 2;
plot (x, y, 'r', "linewidth", 2);
title ("S");
xlabel ("x");
ylabel ("y");

x = linspace (-10, 10, 100);
y = x .^ 2;
plot (x, y, 'r', "linewidth", 2);
title ("S");
xlabel ("x");
ylabel ("y");
axis ([-5, 5, -5, 20]);


How to filter data to remove noise in MATLAB? I have 2 arrays of 800000 input and output data samples. The system is a kind of oven that works from 0 to 10 volts. The sample time is 0.001s.

I have to identify the model of this system, but first of all, given that the data are clearly dirty, I would like to filter the noise.

How can I do it with the System Identification Toolbox of MATLAB?
Moreover, how can I estimate the cutoff frequency to remove the noise?

EDIT:

As suggested, here below are the sampled data plot. As you can see the voltage inputs looks good, but the temperature outputs are affect of noise. I would like to remove it, in order have the output more "smooth".

You can try to remove noise (especially occasional spikes) by non-linear filter. I suggest to use median filter (http://en.wikipedia.org/wiki/Median_filter), as your spikes has length only few points (as I see). There is matlab function - medfilt1. It is probably that such filtration will be enough.


Compensating for measurement errors I have a system where I sample some data periodically (every 10usec). The shape of the data is triangle, in other words it linearly increases and decreases in time. (Both theory and practice are in agreement, I plot the ADC output and I clearly see linear increase and decrease although there is some non-linearity but let's ignore that for the sake of discussion) However the frequency of the signal is not fixed. It may increase very fast or very slow. (narrow or wide triangle)

I am trying to locate the true peak point of the data with 500nsec accuracy, however my ADC samples every 10usec. One item on my favor is that the data values increase and decrease in sequence. I get worst case 50 samples, best case 1000 samples. (in worst case, I will have 25 samples in increasing order and 25 in decreasing order)

How can I approach this issue using DSP techniques? One idea I had was to calculate peak point from rise and fall angles. (I basically extend the left and right side of the triangle and look at where those lines cross each other). I cannot think of any other way. I welcome suggestions.

If you can reliably determine which samples are on a rising portion and which are on a falling portion, then perhaps you can do a linear regression on each (near)monotonic segment, and solve for the intersections of the computed regression lines.

If you know the non-linearity equation, then you can do a regression against that curve instead of using linear regression, and compute or numerically solve for the curve intersections, if needed.

This regression method might also help with error or robustness estimates, if needed.


Confused in Smoothing a signal in the Frequency domain First of all i am new in DSP, so i have no solid education in the field.

I have convert my time domain acceleration-seismic data to frequency using DFT and i am trying to smooth the DFT data, in order to have a more appropriate view of them. I didn't know how exactly to approach my problem until i found a python code in the Scipy-Cookbook that smooths signals with convolution of them with a specific window (e.g. hanning, hamming).

You can find the code here: https://scipy-cookbook.readthedocs.io/items/SignalSmooth.html

I tried this code and the results seem quite promising as you can see.

So my questions are:
1) Is this a theoretically accepted method to smooth data in DFT domain or only applies in time domain ?  (references are welcome!)
2) Are there different methods to smooth DFT data?

Thank you!

1) Is this a theoretically accepted method to smooth data in DFT domain or only applies in time domain ?

Sort of. It certainly can be done this way. Doing it directly in the DFT (instead of the PSD) is risky since it's a complex number you can get a lot of cancellation if the phase is fluctuating a lot.

You also need to decide how you manage the "edges": do you do circular or linear convolution

  2) Are there different methods to smooth DFT data?

Yes. There are as many methods as there are different sets of requirements and application. You plotted your data on a log frequency scale, so you may be better off with "1/Nth octave smoothing algorithm"


Complex Gaussian Channel + AWGN noise If I define a channel to have complex Gaussian coefficients, do I still need to add noise to the signal passing through this channel using the command "awgn" or would that be redundant?

Yes, you still need to add noise to the signal. When you have complex Gaussian coefficients, it says about the probability distribution of the filter coefficients. You still need to add the noise after the filtering process.


Fourier series - time shift and scaling What will be the new Fourier series coefficients when we shift and scale a periodic signal?
Scaling alone will only affect fundamental frequency. But how to calculate new coefficients of shifted and scaled version.
I tried searching, but couldn't find an answer where both properties are used.
Please help.

Consider a continuous-time periodic signal $x_1(t)$ whose fundamental period is $T_1$, fundamental radian frequency is $\omega_1 = \frac{2\pi}{T_1}$ and CTFS (continuous-time Fourier series) coefficients are: $$ x_1 \longleftrightarrow c_k = \frac{1}{T_1} \int_{<T_1>} x_1(t) e^{-j\frac{2\pi k}{T_1}t} dt ~~ , ~~\text{for}~~ k=0,\pm 1, \pm 2...$$

We want to find the new CTFS coefficients, denoted as $d_k$, associated with the new signal $x_2(t)$ , with period $T_2$ and fundamental radian frequency $\omega_2 = \frac{2\pi}{T_2}$, and related to $x_1(t)$ as follows:
$$x_2(t) = x_1( a(t-b) )$$ where $a$ and $b$ are constants, assuming $a > 0$. Then we say ; $$ x_2 \longleftrightarrow d_k $$

We shall proceed from the above definition to show the relation between the coefficients $d_k$ and $c_k$ as follows:
$$ d_k = \frac{1}{T_2} \int_{-T_2/2}^{T_2/2} x_2(t) e^{-j\frac{2\pi}{T_2} k t} dt $$

First note that $T_2 = T_1/a$. And replacing $x_2$ with $x_1$ yields the following:

$$ d_k = \frac{a}{T_1} \int_{t=-T_1/{2a}}^{t=T_1/{2a}} x_1(a(t-b)) e^{-j\frac{2\pi a}{T_1} k t} dt $$
Now apply a substitution $t-b = t'$, and then replace $t'$ with $t$ again (note that $dt'=dt$)

$$ d_k = \frac{a}{T_1} \int_{t=-\frac{T_1}{2a}-b}^{t=\frac{T_1}{2a}-b} x_1(at) e^{-j\frac{2\pi a}{T_1} k (t+b)} dt $$ We know that for a periodic signal, as long as the integration involves a single period, we can redefine its limits:

$$ d_k = e^{-j \frac{2\pi a}{T_1} k b} \left( \frac{a}{T_1} \int_{t=-T_1/{2a}}^{t=T_1/{2a}} x_1(at) e^{-j\frac{2\pi a}{T_1} k t} dt \right)$$
Now apply substitution $at = t'$ (note that $dt' = a dt$) and then replace $t'$ with $t$ again which yields:

$$ d_k = e^{-j \frac{2\pi a}{T_1} k b} \left( \frac{1}{T_1} \int_{t=-T_1/{2}}^{t=T_1/{2}} x_1(t) e^{-j\frac{2\pi}{T_1} k t} dt \right) $$
Finally recognise that the integral equals $T_1 c_k$ and the relationship between $d_k$ and $c_k$ is:

$$ d_k = e^{-j \frac{2\pi}{T_1} a b k} c_k $$

Note that the coefficients $d_k$ are for a periodic signal with fundamental frequency $\omega_2 = \frac{2\pi}{T_2}$, whereas $c_k$ are for the frequency $\omega_1 = \frac{2\pi}{T_1}$. So this means that eventhough, for example, $|d_3|=|c_3|$, i.e., both coefficients have the same magnitude for $k=3$, those coefficients are placed at different frequencies; $d_3$ is associated with $\omega = \frac{2\pi}{T_1}6$ whereas $c_3$ is associated with the frequency $\omega=\frac{2\pi}{T_1}3$... It can be seen that the new coefficients $d_k$ have the same magnitudes, i.e., $|d_k| = |c_k|$, but different phases.


A Sampling theorem for power signals Any function $f \in PW_{\pi}$ (where $PW_{\pi}$ is the Paley-Wiener space), can be expanded in terms of the orthonormal basis $\{e^{i n x}\}_{n\in \mathbb Z}$ as
$$\hat{f}(x)=\sum_{n\in \mathbb Z}\left\langle\hat{f},e^{i n x}\right\rangle_{L^2(-\pi,\pi)} e^{i n x}$$
where $\langle g, h\rangle_{L^2(-\pi,\pi)}=\frac{1}{2\pi}\int_{-\pi}^\pi g(x) \overline{h(x)} dx$. Taking the inverse Fourier transform in previous equation, we obtain the Whittaker-Kotelnikov-Shannon (WKS) sampling theorem,
$$f(x)=\sum_{n\in \mathbb Z}f(n) \operatorname{sinc}(x-n), \ \ \ x\in \mathbb R$$
where $\operatorname{sinc}(x)$ is the normalized sinc function. My question is the following. Sampling theorem, at least in the version that I wrote above, is related to $L^2$ space and so, to scalar product:
$$\langle g, h\rangle_{L^2(-\pi,\pi)}=\frac{1}{2\pi}\int_{-\pi}^\pi g(x) \overline{h(x)} dx$$
Instead, usually, for power signals, we define the inner product to be
$$\left\langle g\,,\,h\right\rangle = \lim_{T\rightarrow \infty} \frac 1 {2T} \int_{-T}^T g(x)\overline{h(x)}\,dx \tag{1}$$
Is there a version of Sampling theorem related to scalar product (1)? Is it used in signal processing?

References and answers are welcome.

(I premise that I'm not an expert on the signal theory, therefore I do apologise if this question isn't much precise.)

Interesting question! The only example I could easily find is this paper.

Casinovi, G., "Sampling and Ergodic Theorems for Weakly Almost Periodic Signals," in Information Theory, IEEE Transactions on, vol.55, no.4, pp.1883-1897, April 2009
doi: 10.1109/TIT.2009.2013021


Effect of DC component on the whole signal - comparison between normalised and non normalised I have a fourier analysis signal as in the picture attached, where red represents the FFT of movement of the hand of a stroke subject and the blue one is the movement of a healthy subject.

I am doing some analysis called Spectral Arc Length, where I will calculate the spectral arc length to compute the smootheness of movement. (check out this paper: On the analysis of movement smoothness by Sivakumar Balasubramanian), where in the metric, the longer the spectral arc length, the less smooth the movement is.

The first image down here is the original data, where we can see that the DC component of the stroke subject is higher than the healthy subject, and the amplitude of the frequency signal is higher, causing the arc length of the signal to be larger, signifies less smooth movement.

MATLAB CODE:

y=variable;
Ts=1/Fs;
L= size(y,2); % Length of signal
NFFT = 2^(ceil(log2(L))+4); % Next power of 2 from length of y
Y = fft( y, NFFT )/L;

f = Fs/2*linspace(0,1,NFFT/2+1);

however, I have a problem with my signal where the metric suggest than we should normalised the signal to the DC component of the signal (Y=Y/max(Y)). By doing this, my signal turn out to be like the second picture, and when I calculate the spectral arc length metric, turns out that stroke 'apparently' have smoother movement than the healthy (due to shorter spectral arc length), which I am pretty sure should't be the case.

My question is, does the normalization to the DC component makes sense? Does the DC component has any effect on the rest of the signal, where larger DC will cause larger amplitude of the signal?
another option for me is to calculate the arc length after the DC signal (starting from the black marker is put on the signal here)
I would also try to do some wavelet analysis if Fourier analysis doesn't work for my calculation...

Non-normalised data

Normalised data Y= Y/max(Y)

As I understand from the description,you are not just using the DC information but also information at other frequencies.

If this is the case:

Generally,Normalizing a signal means subtracting the mean value of the signal from the data.

Y = Y - mean(Y) .......(1)

Fourier transform is a linear transform.So,any constant scaling proportionately scales up the transform.

Y = Y/max(Y)     .... (2)

If (2) is performed then the transform is scaled by max(Y) but the scaling varies with the data and this can be clearly seen in the 2nd plot of your description.